!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

 module ecosys_mod

!BOP
! !MODULE: ecosys_mod
!
! !DESCRIPTION:
!
!  Multispecies ecosystem based on Doney et al. 1996, Moore et al., 2002
!  Based on POP Global NCAR Nitrogen Ecosystem Model
!  version 0.0 (June 15th, 1998) from S.C. Doney.
!  Based on Doney et al., 1996 model.
!  Climate and Global Dynamics, NCAR
!  (doney@whoi.edu)
!
!  Version 1.0
!  Multispecies, multiple limiting nutrient version of ecosystem
!  based on mixed layer model of Moore et al.(2002).  Implemented here with
!  fixed elemental ratios and including only the diatoms and small
!  phytoplankton, with a parameterization of calcification,
!  by Keith Lindsay and Keith Moore, Fall 2001 - Spring 2002.
!  Calcification parameterization based on Moore et al. 2002.
!
!  Version 2.0, January 2003
!    Adds diazotrophs as a phytoplankton group, (based on Moore et al., 2002a)
!    Allows for variable fe/C for all phytoplankton groups
!     Allows for variable si/C for the diatoms
!     Adds explicit tracers for DON, DOP, DOFe
!     variable remin length scale for detrital soft POM and bSi f(temperature)
!     Extensive modifications to iron scavenging parameterization
!     Addition of a sedimentary dissolved iron source,
!        (implemented in ballast code as excess remin in bottom cell)
!        coded by J.K. Moore, (jkmoore@uci.edu)
!
!   Version 2.01. March 2003
!     corrected O2 bug
!     corrected grazing parameter z_grz bug at depth
!     dust dissolution at depth releases iron,
!     increased length scale for dust diss., increased hard fraction dust
!     no deep ocean reduction in scavenging rates,
!     increase bSi OC/ballast ratio 0.3 -> 0.35,
!     corrected bug in diazotroph photoadaptation, and diat and sp adapatation
!
!   Version 2.02.
!     corrected bug in Fe_scavenge (units for dust), May 2003
!     changed C/N/P ratios to 117/16/1 (Anderson & Sarmiento, 1994)
!
!   Version 2.03., July 2003
!     Remin of DOM no longer temperature dependent,
!     new iron scavenging parameterization added,
!     some dissolution of hard fraction of ballast materials added
!
!   Version 2.1, September 2003
!     modfied iron scavenging and dust dissolution at depth
!
!   Version 2.11, March 2004
!     fixed bug in iron scavenging code, replace dust and POC flux_in w/ flux_out
!
!   Version 2.12, April 2004 - Final version for GBC paper revision,
!     (Questions/comments, Keith Moore - jkmoore@uci.edu
!
!   References
!   Doney, S.C., Glover, D.M., Najjar, R.G., 1996. A new coupled, one-dimensional
!   biological-physical model for the upper ocean: applications to the JGOFS
!   Bermuda Time-Series Study (BATS) site. Deep-Sea Res. II, 43: 591-624.
!
!   Moore, JK, Doney, SC, Kleypas, JA, Glover, DM, Fung, IY, 2002. An intermediate
!   complexity marine ecosystem model for the global domain. Deep-Sea Res. II, 49:
!   403-462.
!
!   Moore, JK, Doney, SC, Glover, DM, Fung, IY, 2002. Iron cycling and nutrient
!   limitation patterns in surface waters of the world ocean. Deep-Sea Res. II,
!   49: 463-507.

! !REVISION HISTORY:

!  SVN:$Id:  $

!-----------------------------------------------------------------------
!  variables/subroutines/function used from other modules
!  The following are used extensively in this ecosys, so are used at
!  the module level. The use statements for variables that are only needed
!  locally are located at the module subprogram level.
!-----------------------------------------------------------------------

! !USES:

   use POP_KindsMod
   use POP_ErrorMod
   use POP_CommMod
   use POP_GridHorzMod
   use POP_FieldMod
   use POP_HaloMod

   use kinds_mod
   use constants
   use communicate
   use broadcast
   use global_reductions
   use blocks
   use domain_size
   use domain
   use exit_mod
   use prognostic
   use grid
   use io
   use io_types
   use io_tools
   use tavg
   use timers
   use passive_tracer_tools
   use named_field_mod
   use forcing_tools
   use time_management
   use ecosys_parms
   use registry
   use named_field_mod
   use co2calc
#ifdef CCSMCOUPLED
   use POP_MCT_vars_mod
   use shr_strdata_mod
#endif

! !INPUT PARAMETERS:
!-----------------------------------------------------------------------
!  include ecosystem parameters
!  all variables from this modules have a parm_ prefix
!-----------------------------------------------------------------------

   implicit none
   save
   private

!-----------------------------------------------------------------------
!  public/private declarations
!-----------------------------------------------------------------------

   public :: &
      ecosys_tracer_cnt,            &
      ecosys_init,                  &
      ecosys_tracer_ref_val,        &
      ecosys_set_sflux,             &
      ecosys_tavg_forcing,          &
      ecosys_set_interior,          &
      ecosys_write_restart,         &
      ecosys_qsw_distrb_const

!-----------------------------------------------------------------------
!  module variables required by forcing_passive_tracer
!-----------------------------------------------------------------------

   integer (int_kind), parameter :: &
      ecosys_tracer_cnt = 24

!-----------------------------------------------------------------------
!  flags controlling which portion of code are executed
!  usefull for debugging
!-----------------------------------------------------------------------

  logical (log_kind) :: &
     lsource_sink, &
     lflux_gas_o2, &
     lflux_gas_co2,&
     locmip_k1_k2_bug_fix

  logical (log_kind), dimension(:,:,:), allocatable :: &
     LAND_MASK

!-----------------------------------------------------------------------
!  relative tracer indices
!-----------------------------------------------------------------------

   integer (int_kind), parameter :: &
      po4_ind          =  1,  & ! dissolved inorganic phosphate
      no3_ind          =  2,  & ! dissolved inorganic nitrate
      sio3_ind         =  3,  & ! dissolved inorganic silicate
      nh4_ind          =  4,  & ! dissolved ammonia
      fe_ind           =  5,  & ! dissolved inorganic iron
      o2_ind           =  6,  & ! dissolved oxygen
      dic_ind          =  7,  & ! dissolved inorganic carbon
      alk_ind          =  8,  & ! alkalinity
      doc_ind          =  9,  & ! dissolved organic carbon
      spC_ind          = 10,  & ! small phytoplankton carbon
      spChl_ind        = 11,  & ! small phytoplankton chlorophyll
      spCaCO3_ind      = 12,  & ! small phytoplankton caco3
      diatC_ind        = 13,  & ! diatom carbon
      diatChl_ind      = 14,  & ! diatom chlorophyll
      zooC_ind         = 15,  & ! zooplankton carbon
      spFe_ind         = 16,  & ! small phytoplankton iron
      diatSi_ind       = 17,  & ! diatom silicon
      diatFe_ind       = 18,  & ! diatom iron
      diazC_ind        = 19,  & ! diazotroph carbon
      diazChl_ind      = 20,  & ! diazotroph Chlorophyll
      diazFe_ind       = 21,  & ! diazotroph iron
      don_ind          = 22,  & ! dissolved organic nitrogen
      dofe_ind         = 23,  & ! dissolved organic iron
      dop_ind          = 24     ! dissolved organic phosphorus

!-----------------------------------------------------------------------
!  derived type & parameter for tracer index lookup
!-----------------------------------------------------------------------

   type(ind_name_pair), dimension(ecosys_tracer_cnt) :: &
      ind_name_table

!   type(ind_name_pair), dimension(ecosys_tracer_cnt) :: &
!      ind_name_table = (/ &
!      ind_name_pair(po4_ind,         'PO4'), &
!      ind_name_pair(no3_ind,         'NO3'), &
!      ind_name_pair(sio3_ind,        'SiO3'), &
!      ind_name_pair(nh4_ind,         'NH4'), &
!      ind_name_pair(fe_ind,          'Fe'), &
!      ind_name_pair(o2_ind,          'O2'), &
!      ind_name_pair(dic_ind,         'DIC'), &
!      ind_name_pair(alk_ind,         'ALK'), &
!      ind_name_pair(doc_ind,         'DOC'), &
!      ind_name_pair(spC_ind,         'spC'), &
!      ind_name_pair(spChl_ind,       'spChl'), &
!      ind_name_pair(spCaCO3_ind,     'spCaCO3'), &
!      ind_name_pair(diatC_ind,       'diatC'), &
!      ind_name_pair(diatChl_ind,     'diatChl'), &
!      ind_name_pair(zooC_ind,        'zooC'), &
!      ind_name_pair(spFe_ind,        'spFe'), &
!      ind_name_pair(diatSi_ind,      'diatSi'), &
!      ind_name_pair(diatFe_ind,      'diatFe'), &
!      ind_name_pair(diazC_ind,       'diazC'), &
!      ind_name_pair(diazChl_ind,     'diazChl'), &
!      ind_name_pair(diazFe_ind,      'diazFe'), &
!      ind_name_pair(don_ind,         'DON'), &
!      ind_name_pair(dofe_ind,        'DOFe'), &
!      ind_name_pair(dop_ind,         'DOP') /)
!
!-----------------------------------------------------------------------
!  options for forcing of gas fluxes
!-----------------------------------------------------------------------

   integer (int_kind), parameter :: &
      gas_flux_forcing_iopt_drv   = 1,   &
      gas_flux_forcing_iopt_file  = 2,   &
      atm_co2_iopt_const          = 1,   &
      atm_co2_iopt_drv_prog       = 2,   &
      atm_co2_iopt_drv_diag       = 3

   integer (int_kind) :: &
      gas_flux_forcing_iopt,             &
      atm_co2_iopt

   real (r8)       :: &
      atm_co2_const  ! value of atmospheric co2 (ppm, dry-air, 1 atm)

   character(char_len) :: &
      gas_flux_forcing_file    ! file containing gas flux forcing fields

!-----------------------------------------------------------------------

   real (r8) :: &
      rest_time_inv_surf,  & ! inverse restoring timescale at surface
      rest_time_inv_deep,  & ! inverse restoring timescale at depth
      rest_z0,             & ! shallow end of transition regime
      rest_z1                ! deep end of transition regime

   type(tracer_read) :: &
      po4_rest,            & ! restoring data for PO4
      no3_rest,            & ! restoring data for NO3
      sio3_rest,           & ! restoring data for SiO3
      gas_flux_fice,       & ! ice fraction for gas fluxes
      gas_flux_ws,         & ! wind speed for gas fluxes
      gas_flux_ap,         & ! atmospheric pressure for gas fluxes
      fesedflux_input        ! namelist input for iron_flux

!-----------------------------------------------------------------------
!  module variables related to ph computations
!-----------------------------------------------------------------------

   real (r8), dimension(:,:,:), allocatable, target :: &
      PH_PREV,        & ! computed ph from previous time step
      IRON_PATCH_FLUX   ! localized iron patch flux

   real (r8), dimension(:,:,:), allocatable :: &
      dust_FLUX_IN      ! dust flux not stored in STF since dust is not prognostic

   real (r8), dimension(:,:,:,:), allocatable :: &
      PH_PREV_3D        ! computed pH_3D from previous time step

!-----------------------------------------------------------------------
!  restoring climatologies for nutrients
!-----------------------------------------------------------------------

   logical (log_kind) :: &
      lrest_po4,  & ! restoring on po4 ?
      lrest_no3,  & ! restoring on no3 ?
      lrest_sio3    ! restoring on sio3 ?

   real (r8), dimension(km) :: &
      nutr_rest_time_inv ! inverse restoring time scale for nutrients (1/secs)

   real (r8), dimension(:,:,:,:), allocatable, target :: &
      PO4_CLIM, NO3_CLIM, SiO3_CLIM

   real (r8), dimension(:,:,:,:), allocatable, target :: &
      FESEDFLUX

   character(char_len) :: &
      nutr_rest_file               ! file containing nutrient fields

!maltrud variable restoring
   logical (log_kind) :: &
      lnutr_variable_restore       ! geographically varying nutrient restoring

   character(char_len) :: &
      nutr_variable_rest_file,   & ! file containing variable restoring info
      nutr_variable_rest_file_fmt  ! format of file containing variable restoring info

   real (r8), dimension(:,:,:), allocatable, target :: &
      NUTR_RESTORE_RTAU            ! inverse restoring timescale for variable
                                   ! interior restoring

   integer (int_kind), dimension(:,:,:), allocatable :: &
      NUTR_RESTORE_MAX_LEVEL       ! maximum level for applying variable
                                   ! interior restoring

   real (r8), dimension(:,:,:,:), allocatable :: &
      INTERP_WORK                  ! temp array for interpolate_forcing output

   type(forcing_monthly_every_ts) :: &
      dust_flux,                 & ! surface dust flux
      iron_flux,                 & ! iron component of surface dust flux
      fice_file,                 & ! ice fraction, if read from file
      xkw_file,                  & ! a * wind-speed ** 2, if read from file
      ap_file                      ! atmoshperic pressure, if read from file

   character(char_len) :: &
      ndep_data_type               ! type of ndep forcing

   type(forcing_monthly_every_ts) :: &
      nox_flux_monthly,          & ! surface NOx species flux, added to nitrate pool
      nhy_flux_monthly             ! surface NHy species flux, added to ammonium pool

   integer (int_kind) :: &
      ndep_shr_stream_year_first, & ! first year in stream to use
      ndep_shr_stream_year_last,  & ! last year in stream to use
      ndep_shr_stream_year_align    ! align ndep_shr_stream_year_first with this model year

   integer (int_kind), parameter :: &
      ndep_shr_stream_var_cnt = 2, & ! number of variables in ndep shr_stream
      ndep_shr_stream_no_ind  = 1, & ! index for NO forcing
      ndep_shr_stream_nh_ind  = 2    ! index for NH forcing

   character(char_len) :: &
      ndep_shr_stream_file          ! file containing domain and input data

   real (r8) :: &
      ndep_shr_stream_scale_factor  ! unit conversion factor

#ifdef CCSMCOUPLED
   type(shr_strdata_type) :: ndep_sdat ! input data stream for ndep
#endif

!-----------------------------------------------------------------------
!  derived type for implicit handling of sinking particulate matter
!-----------------------------------------------------------------------

   type sinking_particle
      real (r8) :: &
         diss,        & ! dissolution length for soft subclass
         gamma,       & ! fraction of production -> hard subclass
         mass,        & ! mass of 1e9 base units in g
         rho            ! QA mass ratio of POC to this particle class

      real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
         sflux_in,    & ! incoming flux of soft subclass (base units/cm^2/sec)
         hflux_in,    & ! incoming flux of hard subclass (base units/cm^2/sec)
         prod,        & ! production term (base units/cm^3/sec)
         sflux_out,   & ! outgoing flux of soft subclass (base units/cm^2/sec)
         hflux_out,   & ! outgoing flux of hard subclass (base units/cm^2/sec)
         remin          ! remineralization term (base units/cm^3/sec)
   end type sinking_particle

!-----------------------------------------------------------------------
!  define tavg id for 2d fields related to surface fluxes
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      tavg_ECOSYS_IFRAC, &! tavg id for ice fraction
      tavg_ECOSYS_IFRAC_2,&! tavg id for duplicate of ice fraction
      tavg_ECOSYS_XKW,   &! tavg id for xkw
      tavg_ECOSYS_XKW_2, &! tavg id for duplicate of xkw
      tavg_ECOSYS_ATM_PRESS, &! tavg id for atmospheric pressure
      tavg_PV_O2,        &! tavg id for o2 piston velocity
      tavg_SCHMIDT_O2,   &! tavg id for O2 schmidt number
      tavg_O2SAT,        &! tavg id for O2 saturation
      tavg_CO2STAR,      &! tavg id for co2star
      tavg_DCO2STAR,     &! tavg id for dco2star
      tavg_pCO2SURF,     &! tavg id for surface pco2
      tavg_DpCO2,        &! tavg id for delta pco2
      tavg_DpCO2_2,      &! tavg id for duplicate of delta pco2
      tavg_PV_CO2,       &! tavg id for co2 piston velocity
      tavg_SCHMIDT_CO2,  &! tavg id for co2 schmidt number
      tavg_DIC_GAS_FLUX, &! tavg id for dic flux
      tavg_DIC_GAS_FLUX_2,&! tavg id for duplicate of dic flux
      tavg_O2_GAS_FLUX_2,&! tavg id for duplicate of O2 flux
      tavg_PH,           &! tavg id for surface pH
      tavg_ATM_CO2,      &! tavg id for atmospheric CO2
      tavg_IRON_FLUX,    &! tavg id for dust flux
      tavg_DUST_FLUX,    &! tavg id for dust flux
      tavg_NOx_FLUX,     &! tavg id for nox flux
      tavg_NHy_FLUX       ! tavg id for nhy flux

!-----------------------------------------------------------------------
!  define tavg id for nonstandard 2d fields
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      tavg_O2_ZMIN,      &! tavg id for vertical minimum of O2
      tavg_O2_ZMIN_DEPTH  ! tavg id for depth of vertical minimum of O2

!-----------------------------------------------------------------------
!  define tavg id for nonstandard 3d fields
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      tavg_O2_PRODUCTION,&! tavg id for o2 production
      tavg_O2_CONSUMPTION,&! tavg id for o2 consumption
      tavg_AOU,          &! tavg id for AOU
      tavg_PO4_RESTORE,  &! tavg id for po4 restoring
      tavg_NO3_RESTORE,  &! tavg id for no3 restoring
      tavg_SiO3_RESTORE, &! tavg id for sio3 restoring
      tavg_PAR_avg,      &! tavg id for available radiation avg over mixed layer
      tavg_POC_FLUX_IN,  &! tavg id for poc flux into cell
      tavg_POC_PROD,     &! tavg id for poc production
      tavg_POC_REMIN,    &! tavg id for poc remineralization
      tavg_POC_ACCUM,    &! tavg id for poc accumulation
      tavg_CaCO3_FLUX_IN,&! tavg id for caco3 flux into cell
      tavg_CaCO3_PROD,   &! tavg id for caco3 production
      tavg_CaCO3_REMIN,  &! tavg id for caco3 remineralization
      tavg_SiO2_FLUX_IN, &! tavg id for sio2 flux into cell
      tavg_SiO2_PROD,    &! tavg id for sio2 production
      tavg_SiO2_REMIN,   &! tavg id for sio2 remineralization
      tavg_dust_FLUX_IN, &! tavg id for dust flux into cell
      tavg_dust_REMIN,   &! tavg id for dust remineralization
      tavg_P_iron_FLUX_IN, &! tavg id for p_iron flux into cell
      tavg_P_iron_PROD,    &! tavg id for p_iron production
      tavg_P_iron_REMIN,   &! tavg id for p_iron remineralization
      tavg_graze_sp,     &! tavg id for small phyto grazing
      tavg_graze_diat,   &! tavg id for diatom grazing
      tavg_graze_diaz,   &! tavg id for diaz grazing
      tavg_graze_TOT,    &! tavg id for total grazing
      tavg_sp_loss,      &! tavg id for small phyto loss
      tavg_diat_loss,    &! tavg id for diatom loss
      tavg_diaz_loss,    &! tavg id for diaz loss
      tavg_zoo_loss,     &! tavg id for zooplankton loss
      tavg_sp_agg,       &! tavg id for small phyto aggregate
      tavg_diat_agg       ! tavg id for diatom aggregate

!-----------------------------------------------------------------------
!  define tavg id for MORE nonstandard 3d fields
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      tavg_photoC_sp,            &! tavg id for small phyto C fixation
      tavg_CaCO3_form,           &! tavg id for CaCO3 formation
      tavg_photoC_diat,          &! tavg id for diatom C fixation
      tavg_photoC_diaz,          &! tavg id for diaz C fixation
      tavg_photoC_TOT,           &! tavg id for total C fixation
      tavg_photoC_sp_zint,       &! tavg id for small phyto C fixation vertical integral
      tavg_CaCO3_form_zint,      &! tavg id for CaCO3 formation vertical integral
      tavg_photoC_diat_zint,     &! tavg id for diatom C fixation vertical integral
      tavg_photoC_diaz_zint,     &! tavg id for diaz C fixation vertical integral
      tavg_photoC_TOT_zint,      &! tavg id for total C fixation vertical integral
      tavg_photoC_NO3_sp,        &! tavg id for small phyto C fixation from NO3
      tavg_photoC_NO3_sp_zint,   &! tavg id for small phyto C fixation from NO3 vertical integral
      tavg_photoC_NO3_diat,      &! tavg id for diatom C fixation from NO3
      tavg_photoC_NO3_diat_zint, &! tavg id for diatom C fixation from NO3 vertical integral
      tavg_photoC_NO3_diaz,      &! tavg id for diaz C fixation from NO3
      tavg_photoC_NO3_diaz_zint, &! tavg id for diaz C fixation from NO3 vertical integral
      tavg_photoC_NO3_TOT,       &! tavg id for total C fixation from NO3
      tavg_photoC_NO3_TOT_zint    ! tavg id for total C fixation from NO3 vertical integral

!-----------------------------------------------------------------------
!  define tavg id for MORE nonstandard 3d fields
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      tavg_DOC_prod,       &! tavg id for doc production
      tavg_DOC_remin,      &! tavg id for doc remineralization
      tavg_DON_prod,       &! tavg id for don production
      tavg_DON_remin,      &! tavg id for don remineralization
      tavg_DOFe_prod,      &! tavg id for dofe production
      tavg_DOFe_remin,     &! tavg id for dofe remineralization
      tavg_DOP_prod,       &! tavg id for dop production
      tavg_DOP_remin,      &! tavg id for dop remineralization
      tavg_Fe_scavenge,    &! tavg id for iron scavenging
      tavg_Fe_scavenge_rate,   &! tavg id for iron scavenging rate
      tavg_sp_N_lim,       &! tavg id for small phyto N limitation
      tavg_sp_Fe_lim,      &! tavg id for small phyto Fe limitation
      tavg_sp_PO4_lim,     &! tavg id for small phyto PO4 limitation
      tavg_sp_light_lim,   &! tavg id for small phyto light limitation
      tavg_diat_N_lim,     &! tavg id for diatom N limitation
      tavg_diat_Fe_lim,    &! tavg id for diatom Fe limitation
      tavg_diat_PO4_lim,   &! tavg id for diatom PO4 limitation
      tavg_diat_SiO3_lim,  &! tavg id for diatom SiO3 limitation
      tavg_diat_light_lim, &! tavg id for diatom light limitation
      tavg_diaz_Fe_lim,    &! tavg id for diaz Fe limitation
      tavg_diaz_P_lim,     &! tavg id for diaz PO4 limitation
      tavg_diaz_light_lim, &! tavg id for diaz light limitation
      tavg_diaz_Nfix,      &! tavg id for diaz N fixation
      tavg_bSi_form,       &! tavg id for diatom Si uptake
      tavg_NITRIF,         &! tavg id for nitrification
      tavg_DENITRIF,       &! tavg id for denitrification
      tavg_photoFe_sp,     &! tavg id for small phyto Fe uptake
      tavg_photoFe_diat,   &! tavg id for diatom Fe uptake
      tavg_photoFe_diaz     ! tavg id for diaz Fe uptake

   integer (int_kind) :: &
      tavg_CO3,            &! tavg id for 3D carbonate ion
      tavg_HCO3,           &! tavg id for 3D bicarbonate ion
      tavg_H2CO3,          &! tavg id for 3D carbonic acid
      tavg_pH_3D,          &! tavg id for 3D pH
      tavg_co3_sat_calc,   &! tavg id for co3 concentration at calcite saturation
      tavg_zsatcalc,       &! tavg id for calcite saturation depth
      tavg_co3_sat_arag,   &! tavg id for co3 concentration at aragonite saturation
      tavg_zsatarag         ! tavg id for aragonite saturation depth

!-----------------------------------------------------------------------
!  define array for holding flux-related quantities that need to be time-averaged
!  this is necessary since the forcing routines are called before tavg flags
!-----------------------------------------------------------------------

   real (r8), dimension(:,:,:,:), allocatable :: &
      ECO_SFLUX_TAVG

!-----------------------------------------------------------------------
!  average surface tracer value related variables
!  used as reference value for virtual flux computations
!-----------------------------------------------------------------------

   logical (log_kind), dimension(ecosys_tracer_cnt) :: &
      vflux_flag                ! which tracers get virtual fluxes applied

   integer (int_kind) :: &
      comp_surf_avg_flag        ! time flag id for computing average
                                ! surface tracer values

   real (r8), dimension(ecosys_tracer_cnt) :: &
      surf_avg                  ! average surface tracer values

   logical (log_kind) :: &
      ecosys_qsw_distrb_const

!-----------------------------------------------------------------------
!  iron patch fertilization
!-----------------------------------------------------------------------

   logical (log_kind) :: &
      liron_patch               ! flag for iron patch fertilization

   character(char_len) :: &
      iron_patch_flux_filename  ! file containing name of iron patch file

   integer (int_kind) :: &
      iron_patch_month          !  integer month to add patch flux

!-----------------------------------------------------------------------
!  timers
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      ecosys_shr_strdata_advance_timer,        &
      ecosys_comp_CO3terms_timer,              &
      ecosys_interior_timer,                   &
      ecosys_sflux_timer

!-----------------------------------------------------------------------
!  named field indices
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      totChl_surf_nf_ind = 0,    & ! total chlorophyll in surface layer
      sflux_co2_nf_ind   = 0,    & ! air-sea co2 gas flux
      atm_co2_nf_ind     = 0       ! atmospheric co2

!-----------------------------------------------------------------------

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
      PAR_out           ! photosynthetically available radiation (W/m^2)

!-----------------------------------------------------------------------

   real (r8), parameter :: &
      phlo_surf_init = 7.0_r8, & ! low bound for surface ph for no prev soln
      phhi_surf_init = 9.0_r8, & ! high bound for surface ph for no prev soln
      phlo_3d_init = 6.0_r8,   & ! low bound for subsurface ph for no prev soln
      phhi_3d_init = 9.0_r8,   & ! high bound for subsurface ph for no prev soln
      del_ph = 0.20_r8           ! delta-ph for prev soln

!*****************************************************************************

contains

!*****************************************************************************
!BOP
! !IROUTINE: ecosys_init
! !INTERFACE:

 subroutine ecosys_init(init_ts_file_fmt, read_restart_filename, &
                        tracer_d_module, TRACER_MODULE, tadvect_ctype, &
                        errorCode)

! !DESCRIPTION:
!  Initialize ecosys tracer module. This involves setting metadata, reading
!  the module namelist, setting initial conditions, setting up forcing,
!  and defining additional tavg variables.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   character (*), intent(in) :: &
      init_ts_file_fmt,    & ! format (bin or nc) for input file
      read_restart_filename  ! file name for restart file

! !INPUT/OUTPUT PARAMETERS:

   type (tracer_field), dimension(ecosys_tracer_cnt), intent(inout) :: &
      tracer_d_module   ! descriptors for each tracer

   real (r8), dimension(nx_block,ny_block,km,ecosys_tracer_cnt,3,max_blocks_clinic), &
      intent(inout) :: TRACER_MODULE

! !OUTPUT PARAMETERS:

   character (char_len), dimension(ecosys_tracer_cnt), intent(out) :: &
      tadvect_ctype     ! advection method for ecosys tracers

   integer (POP_i4), intent(out) :: &
      errorCode

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   character(*), parameter :: subname = 'ecosys_mod:ecosys_init'

   character(char_len) :: &
      init_ecosys_option,        & ! option for initialization of bgc
      init_ecosys_init_file,     & ! filename for option 'file'
      init_ecosys_init_file_fmt, & ! file format for option 'file'
      comp_surf_avg_freq_opt,    & ! choice for freq of comp_surf_avg
      gas_flux_forcing_opt,      & ! option for forcing gas fluxes
      atm_co2_opt,               & ! option for atmospheric co2 concentration
      ecosys_tadvect_ctype         ! advection method for ecosys tracers

   type(tracer_read), dimension(ecosys_tracer_cnt) :: &
      tracer_init_ext              ! namelist variable for initializing tracers

   type(tracer_read) :: &
      dust_flux_input,           & ! namelist input for dust_flux
      iron_flux_input,           & ! namelist input for iron_flux
      nox_flux_monthly_input,    & ! namelist input for nox_flux_monthly
      nhy_flux_monthly_input       ! namelist input for nhy_flux_monthly

   logical (log_kind) :: &
      default,                   & ! arg to init_time_flag
      lnml_found,                & ! Was ecosys_nml found ?
      lmarginal_seas               ! Is ecosystem active in marginal seas ?

   integer (int_kind) :: &
      n,                         & ! index for looping over tracers
      k,                         & ! index for looping over depth levels
      l,                         & ! index for looping over time levels
      ind,                       & ! tracer index for tracer name from namelist
      iblock,                    & ! index for looping over blocks
      nml_error                    ! namelist i/o error flag

   integer (int_kind) :: &
      freq_opt, freq,            & ! args for init_time_flag
      comp_surf_avg_freq_iopt,   & ! choice for freq of comp_surf_avg
      comp_surf_avg_freq           ! choice for freq of comp_surf_avg

   logical (log_kind) :: &
      use_nml_surf_vals            ! do namelist surf values override values from restart file

   logical (log_kind) :: &
      lecovars_full_depth_tavg     ! should ecosystem vars be written full depth

!-----------------------------------------------------------------------
!  values to be used when comp_surf_avg_freq_opt==never
!-----------------------------------------------------------------------

   real (r8) :: &
      surf_avg_dic_const, surf_avg_alk_const

   namelist /ecosys_nml/ &
      init_ecosys_option, init_ecosys_init_file, tracer_init_ext, &
      init_ecosys_init_file_fmt, &
      dust_flux_input, iron_flux_input, fesedflux_input, &
      ndep_data_type, nox_flux_monthly_input, nhy_flux_monthly_input, &
      ndep_shr_stream_year_first, ndep_shr_stream_year_last, &
      ndep_shr_stream_year_align, ndep_shr_stream_file, &
      ndep_shr_stream_scale_factor, &
      gas_flux_forcing_opt, gas_flux_forcing_file, &
      gas_flux_fice, gas_flux_ws, gas_flux_ap, &
      lrest_po4, lrest_no3, lrest_sio3, &
      rest_time_inv_surf, rest_time_inv_deep, rest_z0, rest_z1, &
      nutr_rest_file, po4_rest, no3_rest, sio3_rest, &
      comp_surf_avg_freq_opt, comp_surf_avg_freq,  &
      use_nml_surf_vals, surf_avg_dic_const, surf_avg_alk_const, &
      ecosys_qsw_distrb_const, lmarginal_seas, &
      lsource_sink, lflux_gas_o2, lflux_gas_co2, locmip_k1_k2_bug_fix, &
      lnutr_variable_restore, nutr_variable_rest_file,  &
      nutr_variable_rest_file_fmt,atm_co2_opt,atm_co2_const, &
      ecosys_tadvect_ctype, &
      liron_patch,iron_patch_flux_filename,iron_patch_month, &
      lecovars_full_depth_tavg

   character (char_len) :: &
      ecosys_restart_filename  ! modified file name for restart file

   real (r8), dimension (nx_block,ny_block) :: WORK

!-----------------------------------------------------------------------
!  initialize name table 
!-----------------------------------------------------------------------

   errorCode = POP_Success

   ind_name_table( 1) = ind_name_pair(po4_ind,     'PO4')
   ind_name_table( 2) = ind_name_pair(no3_ind,     'NO3')
   ind_name_table( 3) = ind_name_pair(sio3_ind,    'SiO3')
   ind_name_table( 4) = ind_name_pair(nh4_ind,     'NH4')
   ind_name_table( 5) = ind_name_pair(fe_ind,      'Fe')
   ind_name_table( 6) = ind_name_pair(o2_ind,      'O2')
   ind_name_table( 7) = ind_name_pair(dic_ind,     'DIC')
   ind_name_table( 8) = ind_name_pair(alk_ind,     'ALK')
   ind_name_table( 9) = ind_name_pair(doc_ind,     'DOC')
   ind_name_table(10) = ind_name_pair(spC_ind,     'spC')
   ind_name_table(11) = ind_name_pair(spChl_ind,   'spChl')
   ind_name_table(12) = ind_name_pair(spCaCO3_ind, 'spCaCO3')
   ind_name_table(13) = ind_name_pair(diatC_ind,   'diatC')
   ind_name_table(14) = ind_name_pair(diatChl_ind, 'diatChl')
   ind_name_table(15) = ind_name_pair(zooC_ind,    'zooC')
   ind_name_table(16) = ind_name_pair(spFe_ind,    'spFe')
   ind_name_table(17) = ind_name_pair(diatSi_ind,  'diatSi')
   ind_name_table(18) = ind_name_pair(diatFe_ind,  'diatFe')
   ind_name_table(19) = ind_name_pair(diazC_ind,   'diazC')
   ind_name_table(20) = ind_name_pair(diazChl_ind, 'diazChl')
   ind_name_table(21) = ind_name_pair(diazFe_ind,  'diazFe')
   ind_name_table(22) = ind_name_pair(don_ind,     'DON')
   ind_name_table(23) = ind_name_pair(dofe_ind,    'DOFe')
   ind_name_table(24) = ind_name_pair(dop_ind,     'DOP')

!-----------------------------------------------------------------------
!  initialize forcing_monthly_every_ts variables
!-----------------------------------------------------------------------

   call init_forcing_monthly_every_ts(dust_flux)
   call init_forcing_monthly_every_ts(iron_flux)
   call init_forcing_monthly_every_ts(fice_file)
   call init_forcing_monthly_every_ts(xkw_file)
   call init_forcing_monthly_every_ts(ap_file)
   call init_forcing_monthly_every_ts(nox_flux_monthly)
   call init_forcing_monthly_every_ts(nhy_flux_monthly)

!-----------------------------------------------------------------------
!  initialize ecosystem parameters
!-----------------------------------------------------------------------

   call ecosys_parms_init

!-----------------------------------------------------------------------
!  initialize tracer_d values
!-----------------------------------------------------------------------

   do n = 1, ecosys_tracer_cnt
      tracer_d_module(n)%short_name = ind_name_table(n)%name
      if (any(n == (/ spChl_ind, diatChl_ind, diazChl_ind /))) then
         tracer_d_module(n)%units      = 'mg/m^3'
         tracer_d_module(n)%tend_units = 'mg/m^3/s'
         tracer_d_module(n)%flux_units = 'mg/m^3 cm/s'
      else if (n == alk_ind) then
         tracer_d_module(n)%units      = 'meq/m^3'
         tracer_d_module(n)%tend_units = 'meq/m^3/s'
         tracer_d_module(n)%flux_units = 'meq/m^3 cm/s'
      else
         tracer_d_module(n)%units      = 'mmol/m^3'
         tracer_d_module(n)%tend_units = 'mmol/m^3/s'
         tracer_d_module(n)%flux_units = 'mmol/m^3 cm/s'
      endif
   end do
   tracer_d_module(po4_ind    )%long_name='Dissolved Inorganic Phosphate'
   tracer_d_module(no3_ind    )%long_name='Dissolved Inorganic Nitrate'
   tracer_d_module(sio3_ind   )%long_name='Dissolved Inorganic Silicate'
   tracer_d_module(nh4_ind    )%long_name='Dissolved Ammonia'
   tracer_d_module(fe_ind     )%long_name='Dissolved Inorganic Iron'
   tracer_d_module(o2_ind     )%long_name='Dissolved Oxygen'
   tracer_d_module(dic_ind    )%long_name='Dissolved Inorganic Carbon'
   tracer_d_module(alk_ind    )%long_name='Alkalinity'
   tracer_d_module(doc_ind    )%long_name='Dissolved Organic Carbon'
   tracer_d_module(spC_ind    )%long_name='Small Phytoplankton Carbon'
   tracer_d_module(spChl_ind  )%long_name='Small phytoplankton Chlorophyll'
   tracer_d_module(spCaCO3_ind)%long_name='Small Phytoplankton CaCO3'
   tracer_d_module(diatC_ind  )%long_name='Diatom Carbon'
   tracer_d_module(diatChl_ind)%long_name='Diatom Chlorophyll'
   tracer_d_module(zooC_ind   )%long_name='Zooplankton Carbon'
   tracer_d_module(spFe_ind   )%long_name='Small Phytoplankton Iron'
   tracer_d_module(diatSi_ind )%long_name='Diatom Silicon'
   tracer_d_module(diatFe_ind )%long_name='Diatom Iron'
   tracer_d_module(diazC_ind  )%long_name='Diazotroph Carbon'
   tracer_d_module(diazChl_ind)%long_name='Diazotroph Chlorophyll'
   tracer_d_module(diazFe_ind )%long_name='Diazotroph Iron'
   tracer_d_module(don_ind    )%long_name='Dissolved Organic Nitrogen'
   tracer_d_module(dofe_ind   )%long_name='Dissolved Organic Iron'
   tracer_d_module(dop_ind    )%long_name='Dissolved Organic Phosphorus'

!-----------------------------------------------------------------------
!  default namelist settings
!-----------------------------------------------------------------------

   init_ecosys_option = 'unknown'
   init_ecosys_init_file = 'unknown'
   init_ecosys_init_file_fmt = 'bin'

   gas_flux_forcing_opt  = 'drv'
   gas_flux_forcing_file = 'unknown'

   gas_flux_fice%filename     = 'unknown'
   gas_flux_fice%file_varname = 'FICE'
   gas_flux_fice%scale_factor = c1
   gas_flux_fice%default_val  = c0
   gas_flux_fice%file_fmt     = 'bin'

   gas_flux_ws%filename     = 'unknown'
   gas_flux_ws%file_varname = 'XKW'
   gas_flux_ws%scale_factor = c1
   gas_flux_ws%default_val  = c0
   gas_flux_ws%file_fmt     = 'bin'

   gas_flux_ap%filename     = 'unknown'
   gas_flux_ap%file_varname = 'P'
   gas_flux_ap%scale_factor = c1
   gas_flux_ap%default_val  = c0
   gas_flux_ap%file_fmt     = 'bin'

   lrest_po4      = .false.
   lrest_no3      = .false.
   lrest_sio3     = .false.
   nutr_rest_file = 'unknown'

   rest_time_inv_surf = c0
   rest_time_inv_deep = c0
   rest_z0            = c1000
   rest_z1            = c2 * c1000

   po4_rest%filename     = 'unknown'
   po4_rest%file_varname = tracer_d_module(po4_ind)%short_name
   po4_rest%scale_factor = c1
   po4_rest%default_val  = c0
   po4_rest%file_fmt     = 'bin'

   no3_rest%filename     = 'unknown'
   no3_rest%file_varname = tracer_d_module(no3_ind)%short_name
   no3_rest%scale_factor = c1
   no3_rest%default_val  = c0
   no3_rest%file_fmt     = 'bin'

   sio3_rest%filename     = 'unknown'
   sio3_rest%file_varname = tracer_d_module(sio3_ind)%short_name
   sio3_rest%scale_factor = c1
   sio3_rest%default_val  = c0
   sio3_rest%file_fmt     = 'bin'

!maltrud variable restoring
   lnutr_variable_restore      = .false.
   nutr_variable_rest_file     = 'unknown'
   nutr_variable_rest_file_fmt = 'bin'

   dust_flux_input%filename     = 'unknown'
   dust_flux_input%file_varname = 'dust_flux'
   dust_flux_input%scale_factor = c1
   dust_flux_input%default_val  = c0
   dust_flux_input%file_fmt     = 'bin'

   iron_flux_input%filename     = 'unknown'
   iron_flux_input%file_varname = 'iron_flux'
   iron_flux_input%scale_factor = c1
   iron_flux_input%default_val  = c0
   iron_flux_input%file_fmt     = 'bin'

   fesedflux_input%filename     = 'unknown'
   fesedflux_input%file_varname = 'FESEDFLUXIN'
   fesedflux_input%scale_factor = c1
   fesedflux_input%default_val  = c0
   fesedflux_input%file_fmt     = 'bin'

   ndep_data_type              = 'monthly-calendar'

   nox_flux_monthly_input%filename     = 'unknown'
   nox_flux_monthly_input%file_varname = 'nox_flux'
   nox_flux_monthly_input%scale_factor = c1
   nox_flux_monthly_input%default_val  = c0
   nox_flux_monthly_input%file_fmt     = 'bin'

   nhy_flux_monthly_input%filename     = 'unknown'
   nhy_flux_monthly_input%file_varname = 'nhy_flux'
   nhy_flux_monthly_input%scale_factor = c1
   nhy_flux_monthly_input%default_val  = c0
   nhy_flux_monthly_input%file_fmt     = 'bin'

   ndep_shr_stream_year_first = 1
   ndep_shr_stream_year_last  = 1
   ndep_shr_stream_year_align = 1
   ndep_shr_stream_file       = 'unknown'
   ndep_shr_stream_scale_factor = c1

   do n = 1,ecosys_tracer_cnt
      tracer_init_ext(n)%mod_varname  = 'unknown'
      tracer_init_ext(n)%filename     = 'unknown'
      tracer_init_ext(n)%file_varname = 'unknown'
      tracer_init_ext(n)%scale_factor = c1
      tracer_init_ext(n)%default_val  = c0
      tracer_init_ext(n)%file_fmt     = 'bin'
   end do

   lmarginal_seas        = .true.
   lsource_sink          = .true.
   lflux_gas_o2          = .true.
   lflux_gas_co2         = .true.
   locmip_k1_k2_bug_fix  = .true.

   comp_surf_avg_freq_opt        = 'never'
   comp_surf_avg_freq            = 1
   use_nml_surf_vals             = .false.
   surf_avg_dic_const            = 1944.0_r8
   surf_avg_alk_const            = 2225.0_r8

   ecosys_qsw_distrb_const  = .true.

   liron_patch              = .false.
   iron_patch_flux_filename = 'unknown_iron_patch_filename'
   iron_patch_month         = 1

   atm_co2_opt   = 'const'
   atm_co2_const = 280.0_r8

   ecosys_tadvect_ctype = 'base_model'

   lecovars_full_depth_tavg = .false.

   if (my_task == master_task) then
      open (nml_in, file=nml_filename, status='old',iostat=nml_error)
      if (nml_error /= 0) then
         nml_error = -1
      else
         nml_error =  1
      endif
      do while (nml_error > 0)
         read(nml_in, nml=ecosys_nml,iostat=nml_error)
      end do
      if (nml_error == 0) close(nml_in)
   endif

   call broadcast_scalar(nml_error, master_task)
   if (nml_error /= 0) then
      call document(subname, 'ecosys_nml not found')
      call exit_POP(sigAbort, 'ERROR : stopping in '/&
                           &/ subname)
   endif

   if (my_task == master_task) then
      write(stdout,blank_fmt)
      write(stdout,ndelim_fmt)
      write(stdout,blank_fmt)
      write(stdout,*) ' ecosys:'
      write(stdout,blank_fmt)
      write(stdout,*) ' ecosys_nml namelist settings:'
      write(stdout,blank_fmt)
      write(stdout,ecosys_nml)
      write(stdout,blank_fmt)
      write(stdout,delim_fmt)
   endif

!-----------------------------------------------------------------------
!  broadcast all namelist variables
!-----------------------------------------------------------------------

   call broadcast_scalar(init_ecosys_option, master_task)
   call broadcast_scalar(init_ecosys_init_file, master_task)
   call broadcast_scalar(init_ecosys_init_file_fmt, master_task)

   call broadcast_scalar(gas_flux_forcing_opt, master_task)
   if (trim(gas_flux_forcing_opt) == 'drv') then
      gas_flux_forcing_iopt = gas_flux_forcing_iopt_drv
   else if (trim(gas_flux_forcing_opt) == 'file') then
      gas_flux_forcing_iopt = gas_flux_forcing_iopt_file
   else
      call document(subname, 'gas_flux_forcing_opt', gas_flux_forcing_opt)
      call exit_POP(sigAbort, 'unknown gas_flux_forcing_opt')
   endif

   call broadcast_scalar(gas_flux_forcing_file, master_task)

   call broadcast_scalar(gas_flux_fice%filename, master_task)
   call broadcast_scalar(gas_flux_fice%file_varname, master_task)
   call broadcast_scalar(gas_flux_fice%scale_factor, master_task)
   call broadcast_scalar(gas_flux_fice%default_val, master_task)
   call broadcast_scalar(gas_flux_fice%file_fmt, master_task)

   fice_file%input = gas_flux_fice

   call broadcast_scalar(gas_flux_ws%filename, master_task)
   call broadcast_scalar(gas_flux_ws%file_varname, master_task)
   call broadcast_scalar(gas_flux_ws%scale_factor, master_task)
   call broadcast_scalar(gas_flux_ws%default_val, master_task)
   call broadcast_scalar(gas_flux_ws%file_fmt, master_task)

   xkw_file%input = gas_flux_ws

   call broadcast_scalar(gas_flux_ap%filename, master_task)
   call broadcast_scalar(gas_flux_ap%file_varname, master_task)
   call broadcast_scalar(gas_flux_ap%scale_factor, master_task)
   call broadcast_scalar(gas_flux_ap%default_val, master_task)
   call broadcast_scalar(gas_flux_ap%file_fmt, master_task)

   ap_file%input = gas_flux_ap

   call broadcast_scalar(lrest_po4, master_task)
   call broadcast_scalar(lrest_no3, master_task)
   call broadcast_scalar(lrest_sio3, master_task)
   call broadcast_scalar(nutr_rest_file, master_task)

   call broadcast_scalar(rest_time_inv_surf, master_task)
   call broadcast_scalar(rest_time_inv_deep, master_task)
   call broadcast_scalar(rest_z0, master_task)
   call broadcast_scalar(rest_z1, master_task)

   call broadcast_scalar(po4_rest%filename, master_task)
   call broadcast_scalar(po4_rest%file_varname, master_task)
   call broadcast_scalar(po4_rest%scale_factor, master_task)
   call broadcast_scalar(po4_rest%default_val, master_task)
   call broadcast_scalar(po4_rest%file_fmt, master_task)

   call broadcast_scalar(no3_rest%filename, master_task)
   call broadcast_scalar(no3_rest%file_varname, master_task)
   call broadcast_scalar(no3_rest%scale_factor, master_task)
   call broadcast_scalar(no3_rest%default_val, master_task)
   call broadcast_scalar(no3_rest%file_fmt, master_task)

   call broadcast_scalar(sio3_rest%filename, master_task)
   call broadcast_scalar(sio3_rest%file_varname, master_task)
   call broadcast_scalar(sio3_rest%scale_factor, master_task)
   call broadcast_scalar(sio3_rest%default_val, master_task)
   call broadcast_scalar(sio3_rest%file_fmt, master_task)

!maltrud variable restoring
   call broadcast_scalar(lnutr_variable_restore, master_task)
   call broadcast_scalar(nutr_variable_rest_file, master_task)
   call broadcast_scalar(nutr_variable_rest_file_fmt, master_task)

   call broadcast_scalar(dust_flux_input%filename, master_task)
   call broadcast_scalar(dust_flux_input%file_varname, master_task)
   call broadcast_scalar(dust_flux_input%scale_factor, master_task)
   call broadcast_scalar(dust_flux_input%default_val, master_task)
   call broadcast_scalar(dust_flux_input%file_fmt, master_task)

   dust_flux%input = dust_flux_input

   call broadcast_scalar(iron_flux_input%filename, master_task)
   call broadcast_scalar(iron_flux_input%file_varname, master_task)
   call broadcast_scalar(iron_flux_input%scale_factor, master_task)
   call broadcast_scalar(iron_flux_input%default_val, master_task)
   call broadcast_scalar(iron_flux_input%file_fmt, master_task)

   iron_flux%input = iron_flux_input

   call broadcast_scalar(fesedflux_input%filename, master_task)
   call broadcast_scalar(fesedflux_input%file_varname, master_task)
   call broadcast_scalar(fesedflux_input%scale_factor, master_task)
   call broadcast_scalar(fesedflux_input%default_val, master_task)
   call broadcast_scalar(fesedflux_input%file_fmt, master_task)

   call broadcast_scalar(ndep_data_type, master_task)

   call broadcast_scalar(nox_flux_monthly_input%filename, master_task)
   call broadcast_scalar(nox_flux_monthly_input%file_varname, master_task)
   call broadcast_scalar(nox_flux_monthly_input%scale_factor, master_task)
   call broadcast_scalar(nox_flux_monthly_input%default_val, master_task)
   call broadcast_scalar(nox_flux_monthly_input%file_fmt, master_task)

   nox_flux_monthly%input = nox_flux_monthly_input

   call broadcast_scalar(nhy_flux_monthly_input%filename, master_task)
   call broadcast_scalar(nhy_flux_monthly_input%file_varname, master_task)
   call broadcast_scalar(nhy_flux_monthly_input%scale_factor, master_task)
   call broadcast_scalar(nhy_flux_monthly_input%default_val, master_task)
   call broadcast_scalar(nhy_flux_monthly_input%file_fmt, master_task)

   nhy_flux_monthly%input = nhy_flux_monthly_input

   call broadcast_scalar(ndep_shr_stream_year_first, master_task)
   call broadcast_scalar(ndep_shr_stream_year_last, master_task)
   call broadcast_scalar(ndep_shr_stream_year_align, master_task)
   call broadcast_scalar(ndep_shr_stream_file, master_task)
   call broadcast_scalar(ndep_shr_stream_scale_factor, master_task)

   do n = 1,ecosys_tracer_cnt
      call broadcast_scalar(tracer_init_ext(n)%mod_varname, master_task)
      call broadcast_scalar(tracer_init_ext(n)%filename, master_task)
      call broadcast_scalar(tracer_init_ext(n)%file_varname, master_task)
      call broadcast_scalar(tracer_init_ext(n)%scale_factor, master_task)
      call broadcast_scalar(tracer_init_ext(n)%default_val, master_task)
      call broadcast_scalar(tracer_init_ext(n)%file_fmt, master_task)
   end do

   call broadcast_scalar(comp_surf_avg_freq_opt, master_task)
   call broadcast_scalar(comp_surf_avg_freq, master_task)
   call broadcast_scalar(use_nml_surf_vals, master_task)
   call broadcast_scalar(surf_avg_dic_const, master_task)
   call broadcast_scalar(surf_avg_alk_const, master_task)

   call broadcast_scalar(ecosys_qsw_distrb_const, master_task)

   call broadcast_scalar(lmarginal_seas, master_task)
   call broadcast_scalar(lsource_sink, master_task)
   call broadcast_scalar(lflux_gas_o2, master_task)
   call broadcast_scalar(lflux_gas_co2, master_task)
   call broadcast_scalar(locmip_k1_k2_bug_fix, master_task)

   call broadcast_scalar(liron_patch, master_task)
   call broadcast_scalar(iron_patch_flux_filename, master_task)
   call broadcast_scalar(iron_patch_month, master_task)

   call broadcast_scalar(atm_co2_opt, master_task)
   call broadcast_scalar(atm_co2_const, master_task)

   call broadcast_scalar(ecosys_tadvect_ctype, master_task)
   tadvect_ctype = ecosys_tadvect_ctype

   call broadcast_scalar(lecovars_full_depth_tavg, master_task)

!-----------------------------------------------------------------------
!  set variables immediately dependent on namelist variables
!-----------------------------------------------------------------------

   select case (comp_surf_avg_freq_opt)
   case ('never')
      comp_surf_avg_freq_iopt = freq_opt_never
   case ('nyear')
      comp_surf_avg_freq_iopt = freq_opt_nyear
   case ('nmonth')
      comp_surf_avg_freq_iopt = freq_opt_nmonth
   case default
      call document(subname, 'comp_surf_avg_freq_opt', comp_surf_avg_freq_opt)
      call exit_POP(sigAbort, 'unknown comp_surf_avg_freq_opt')
   end select

  call init_time_flag('ecosys_comp_surf_avg', comp_surf_avg_flag, &
     default=.false., freq_opt=comp_surf_avg_freq_iopt,  &
     freq=comp_surf_avg_freq, owner='ecosys_init')

   select case (atm_co2_opt)
   case ('const')
      atm_co2_iopt = atm_co2_iopt_const
   case ('drv_prog')
      atm_co2_iopt = atm_co2_iopt_drv_prog
   case ('drv_diag')
      atm_co2_iopt = atm_co2_iopt_drv_diag
   case default
      call document(subname, 'atm_co2_opt', atm_co2_opt)
      call exit_POP(sigAbort, 'unknown atm_co2_opt')
   end select

!-----------------------------------------------------------------------
!  namelist consistency checking
!-----------------------------------------------------------------------

   if (use_nml_surf_vals .and. comp_surf_avg_freq_iopt /= freq_opt_never) then
      call document(subname, 'use_nml_surf_vals', use_nml_surf_vals)
      call document(subname, 'comp_surf_avg_freq_opt', comp_surf_avg_freq_opt)
      call exit_POP(sigAbort, 'use_nml_surf_vals can only be .true. if ' /&
                           &/ ' comp_surf_avg_freq_opt is never')
   endif

!-----------------------------------------------------------------------
!  initialize virtual flux flag array
!-----------------------------------------------------------------------

   vflux_flag = .false.
   vflux_flag(dic_ind) = .true.
   vflux_flag(alk_ind) = .true.

!-----------------------------------------------------------------------
!  allocate various ecosys allocatable module variables
!-----------------------------------------------------------------------

   allocate( PH_PREV(nx_block,ny_block,max_blocks_clinic) )
   allocate( dust_FLUX_IN(nx_block,ny_block,max_blocks_clinic) )
   allocate( PH_PREV_3D(nx_block,ny_block,km,max_blocks_clinic) )
   PH_PREV_3D = c0

!-----------------------------------------------------------------------
!  allocate and initialize LAND_MASK
!-----------------------------------------------------------------------

   allocate( LAND_MASK(nx_block,ny_block,nblocks_clinic) )

   if (lmarginal_seas) then
      LAND_MASK = REGION_MASK /= 0
   else
      LAND_MASK = REGION_MASK > 0
   endif

!-----------------------------------------------------------------------
!  initialize tracers
!-----------------------------------------------------------------------

   select case (init_ecosys_option)

   case ('restart', 'ccsm_continue', 'ccsm_branch', 'ccsm_hybrid')

      ecosys_restart_filename = char_blank

      if (init_ecosys_init_file == 'same_as_TS') then
         if (read_restart_filename == 'undefined') then
            call document(subname, 'no restart file to read ecosys from')
            call exit_POP(sigAbort, 'stopping in ' /&
                                 &/ subname)
         endif
         ecosys_restart_filename = read_restart_filename
         init_ecosys_init_file_fmt = init_ts_file_fmt

      else  ! do not read from TS restart file

         ecosys_restart_filename = trim(init_ecosys_init_file)

      endif

      call rest_read_tracer_block(init_ecosys_init_file_fmt, &
                                  ecosys_restart_filename,   &
                                  tracer_d_module,           &
                                  TRACER_MODULE)

      call read_field(init_ecosys_init_file_fmt, &
                      ecosys_restart_filename,   &
                      'PH_SURF', PH_PREV)

      if (use_nml_surf_vals) then
         surf_avg = c0
         surf_avg(dic_ind) = surf_avg_dic_const
         surf_avg(alk_ind) = surf_avg_alk_const
      else
         call extract_surf_avg(init_ecosys_init_file_fmt, &
                               ecosys_restart_filename)
      endif

      call eval_time_flag(comp_surf_avg_flag) ! evaluates time_flag(comp_surf_avg_flag)%value via time_to_do

      if (check_time_flag(comp_surf_avg_flag)) &
         call comp_surf_avg(TRACER_MODULE(:,:,1,:,oldtime,:), &
                            TRACER_MODULE(:,:,1,:,curtime,:))

   case ('file', 'ccsm_startup')
      call document(subname, 'ecosystem vars being read from separate files')

      call file_read_tracer_block(init_ecosys_init_file_fmt, &
                                  init_ecosys_init_file,     &
                                  tracer_d_module,           &
                                  ind_name_table,            &
                                  tracer_init_ext,           &
                                  TRACER_MODULE)

      if (n_topo_smooth > 0) then
         do n = 1, ecosys_tracer_cnt
            do k=1,km
               call fill_points(k,TRACER_MODULE(:,:,k,n,oldtime,:), &
                                errorCode)

               if (errorCode /= POP_Success) then
                  call POP_ErrorSet(errorCode, &
                     'ecosys_init: error in fill points for tracers(oldtime)')
                  return
               endif

               call fill_points(k,TRACER_MODULE(:,:,k,n,curtime,:), &
                                errorCode)

               if (errorCode /= POP_Success) then
                  call POP_ErrorSet(errorCode, &
                     'ecosys_init: error in fill points for tracers(newtime)')
                  return
               endif

            enddo
         enddo
      endif

      PH_PREV = c0

      if (use_nml_surf_vals) then
         surf_avg = c0
         surf_avg(dic_ind) = surf_avg_dic_const
         surf_avg(alk_ind) = surf_avg_alk_const
      else
         call comp_surf_avg(TRACER_MODULE(:,:,1,:,oldtime,:), &
                            TRACER_MODULE(:,:,1,:,curtime,:))
      endif

   case default
      call document(subname, 'init_ecosys_option', init_ecosys_option)
      call exit_POP(sigAbort, 'unknown init_ecosys_option')

   end select

!-----------------------------------------------------------------------
!  register Chl field for short-wave absorption
!  apply land mask to tracers
!  set Chl field for short-wave absorption
!-----------------------------------------------------------------------

   call named_field_register('model_chlorophyll', totChl_surf_nf_ind)

   !$OMP PARALLEL DO PRIVATE(iblock,n,k,WORK)
   do iblock=1,nblocks_clinic
      do n = 1,ecosys_tracer_cnt
         do k = 1,km
            where (.not. LAND_MASK(:,:,iblock) .or. k > KMT(:,:,iblock))
               TRACER_MODULE(:,:,k,n,curtime,iblock) = c0
               TRACER_MODULE(:,:,k,n,oldtime,iblock) = c0
            end where
         end do
      end do

      WORK = max(c0,p5*(TRACER_MODULE(:,:,1,spChl_ind,oldtime,iblock) + &
                        TRACER_MODULE(:,:,1,spChl_ind,curtime,iblock))) + &
             max(c0,p5*(TRACER_MODULE(:,:,1,diatChl_ind,oldtime,iblock) + &
                        TRACER_MODULE(:,:,1,diatChl_ind,curtime,iblock))) + &
             max(c0,p5*(TRACER_MODULE(:,:,1,diazChl_ind,oldtime,iblock) + &
                        TRACER_MODULE(:,:,1,diazChl_ind,curtime,iblock)))
      call named_field_set(totChl_surf_nf_ind, iblock, WORK)
   enddo
   !$OMP END PARALLEL DO

!-----------------------------------------------------------------------
!  timer init
!-----------------------------------------------------------------------

   call get_timer(ecosys_comp_CO3terms_timer, 'comp_CO3terms', &
                  nblocks_clinic, distrb_clinic%nprocs)
   call get_timer(ecosys_interior_timer, 'ECOSYS_INTERIOR', &
                  nblocks_clinic, distrb_clinic%nprocs)
   call get_timer(ecosys_sflux_timer, 'ECOSYS_SFLUX',1, &
                  distrb_clinic%nprocs)
   if (ndep_data_type == 'shr_stream') then
      call get_timer(ecosys_shr_strdata_advance_timer, &
                     'ecosys_shr_strdata_advance',1, distrb_clinic%nprocs)
   endif

!-----------------------------------------------------------------------
!  call other initialization subroutines
!-----------------------------------------------------------------------

   call ecosys_init_tavg
   call ecosys_init_sflux
   call ecosys_init_interior_restore

!-----------------------------------------------------------------------
!  set lfull_depth_tavg flag for short-lived ecosystem tracers
!-----------------------------------------------------------------------

   tracer_d_module(spC_ind    )%lfull_depth_tavg = lecovars_full_depth_tavg
   tracer_d_module(spChl_ind  )%lfull_depth_tavg = lecovars_full_depth_tavg
   tracer_d_module(spCaCO3_ind)%lfull_depth_tavg = lecovars_full_depth_tavg
   tracer_d_module(diatC_ind  )%lfull_depth_tavg = lecovars_full_depth_tavg
   tracer_d_module(diatChl_ind)%lfull_depth_tavg = lecovars_full_depth_tavg
   tracer_d_module(zooC_ind   )%lfull_depth_tavg = lecovars_full_depth_tavg
   tracer_d_module(spFe_ind   )%lfull_depth_tavg = lecovars_full_depth_tavg
   tracer_d_module(diatSi_ind )%lfull_depth_tavg = lecovars_full_depth_tavg
   tracer_d_module(diatFe_ind )%lfull_depth_tavg = lecovars_full_depth_tavg
   tracer_d_module(diazC_ind  )%lfull_depth_tavg = lecovars_full_depth_tavg
   tracer_d_module(diazChl_ind)%lfull_depth_tavg = lecovars_full_depth_tavg
   tracer_d_module(diazFe_ind )%lfull_depth_tavg = lecovars_full_depth_tavg

!-----------------------------------------------------------------------
!EOC

 end subroutine ecosys_init

!***********************************************************************
!BOP
! !IROUTINE: extract_surf_avg
! !INTERFACE:

 subroutine extract_surf_avg(init_ecosys_init_file_fmt, &
                             ecosys_restart_filename)

! !DESCRIPTION:
!  Extract average surface values from restart file.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   character (*), intent(in) :: &
      init_ecosys_init_file_fmt, & ! file format (bin or nc)
      ecosys_restart_filename      ! file name for restart file

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   type (datafile) ::&
      restart_file    ! io file descriptor

   integer (int_kind) :: &
      n               ! tracer index

   character (char_len) :: &
      short_name      ! tracer name temporaries

!-----------------------------------------------------------------------

   surf_avg = c0

   restart_file = construct_file(init_ecosys_init_file_fmt, &
                                 full_name=trim(ecosys_restart_filename), &
                                 record_length=rec_type_dbl, &
                                 recl_words=nx_global*ny_global)

   do n = 1, ecosys_tracer_cnt
      if (vflux_flag(n)) then
         short_name = 'surf_avg_' /&
                   &/ ind_name_table(n)%name
         call add_attrib_file(restart_file, trim(short_name), surf_avg(n))
      endif
   end do

   call data_set(restart_file, 'open_read')

   do n = 1, ecosys_tracer_cnt
      if (vflux_flag(n)) then
         short_name = 'surf_avg_' /&
                   &/ ind_name_table(n)%name
         call extract_attrib_file(restart_file, trim(short_name), surf_avg(n))
      endif
   end do

   call data_set (restart_file, 'close')

   call destroy_file (restart_file)

!-----------------------------------------------------------------------
!EOC

 end subroutine extract_surf_avg

!***********************************************************************
!BOP
! !IROUTINE: ecosys_init_tavg
! !INTERFACE:

 subroutine ecosys_init_tavg

! !DESCRIPTION:
!  call define_tavg_field for nonstandard tavg fields
!
! !REVISION HISTORY:
!  same as module
!

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      var_cnt             ! how many tavg variables are defined

!-----------------------------------------------------------------------
!  2D fields related to surface fluxes
!-----------------------------------------------------------------------

   var_cnt = 0

   call define_tavg_field(tavg_ECOSYS_IFRAC,'ECOSYS_IFRAC',2,          &
                          long_name='Ice Fraction for ecosys fluxes',  &
                          units='fraction', grid_loc='2110',           &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_ECOSYS_IFRAC_2,'ECOSYS_IFRAC_2',2,      &
                          long_name='Ice Fraction for ecosys fluxes',  &
                          units='fraction', grid_loc='2110',           &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_ECOSYS_XKW,'ECOSYS_XKW',2,              &
                          long_name='XKW for ecosys fluxes',           &
                          units='cm/s', grid_loc='2110',               &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_ECOSYS_XKW_2,'ECOSYS_XKW_2',2,          &
                          long_name='XKW for ecosys fluxes',           &
                          units='cm/s', grid_loc='2110',               &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_ECOSYS_ATM_PRESS,'ECOSYS_ATM_PRESS',2,  &
                          long_name='Atmospheric Pressure for ecosys fluxes', &
                          units='atmospheres', grid_loc='2110',        &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_PV_O2,'PV_O2',2,                        &
                          long_name='PV_O2',                           &
                          units='cm/s', grid_loc='2110',               &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_SCHMIDT_O2,'SCHMIDT_O2',2,              &
                          long_name='O2 Schmidt Number',               &
                          units='none', grid_loc='2110',               &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_O2SAT,'O2SAT',2,                        &
                          long_name='O2 Saturation',                   &
                          units='mmol/m^3', grid_loc='2110',           &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_CO2STAR,'CO2STAR',2,                    &
                          long_name='CO2 Star',                        &
                          units='mmol/m^3', grid_loc='2110',           &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_DCO2STAR,'DCO2STAR',2,                  &
                          long_name='D CO2 Star',                      &
                          units='mmol/m^3', grid_loc='2110',           &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_pCO2SURF,'pCO2SURF',2,                  &
                          long_name='surface pCO2',                    &
                          units='ppmv', grid_loc='2110',               &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_DpCO2,'DpCO2',2,                        &
                          long_name='D pCO2',                          &
                          units='ppmv', grid_loc='2110',               &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_DpCO2_2,'DpCO2_2',2,                    &
                          long_name='D pCO2',                          &
                          units='ppmv', grid_loc='2110',               &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_PV_CO2,'PV_CO2',2,                      &
                          long_name='CO2 Piston Velocity',             &
                          units='cm/s', grid_loc='2110',               &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_SCHMIDT_CO2,'SCHMIDT_CO2',2,            &
                          long_name='CO2 Schmidt Number',              &
                          units='none', grid_loc='2110',               &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_DIC_GAS_FLUX,'FG_CO2',2,                &
                          long_name='DIC Surface Gas Flux',            &
                          units='mmol/m^3 cm/s', grid_loc='2110',      &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_DIC_GAS_FLUX_2,'FG_CO2_2',2,            &
                          long_name='DIC Surface Gas Flux',            &
                          units='mmol/m^3 cm/s', grid_loc='2110',      &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_PH,'PH',2,                              &
                          long_name='Surface pH',                      &
                          units='none', grid_loc='2110',               &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_ATM_CO2,'ATM_CO2',2,                    &
                          long_name='Atmospheric CO2',                 &
                          units='ppmv', grid_loc='2110',               &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_O2_GAS_FLUX_2,'STF_O2_2',2,             &
                          long_name='Dissolved Oxygen Surface Flux',   &
                          units='mmol/m^3 cm/s', grid_loc='2110',      &
                          coordinates='TLONG TLAT time')
   var_cnt = var_cnt+1

!-----------------------------------------------------------------------
 
   allocate(ECO_SFLUX_TAVG(nx_block,ny_block,var_cnt,max_blocks_clinic))
   ECO_SFLUX_TAVG = c0
 
!-----------------------------------------------------------------------
!  The follow surface flux related fields are not with the above because
!  their values are directly available from STF in ecosys_tavg_forcing
!  and do not need to be stored in ECO_SFLUX_TAVG.
!-----------------------------------------------------------------------

   call define_tavg_field(tavg_IRON_FLUX,'IRON_FLUX',2,                &
                          long_name='Iron Flux',                       &
                          units='mmol/m^2/s', grid_loc='2110',        &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_DUST_FLUX,'DUST_FLUX',2,                &
                          long_name='Dust Flux',                       &
                          units='g/cm^2/s', grid_loc='2110',           &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_NOx_FLUX,'NOx_FLUX',2,                  &
                          long_name='Flux of NOx from Atmosphere',     &
                          units='nmol/cm^2/s', grid_loc='2110',        &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_NHy_FLUX,'NHy_FLUX',2,                  &
                          long_name='Flux of NHy from Atmosphere',     &
                          units='nmol/cm^2/s', grid_loc='2110',        &
                          coordinates='TLONG TLAT time')

!-----------------------------------------------------------------------
!  nonstandard 2D fields
!-----------------------------------------------------------------------

   call define_tavg_field(tavg_O2_ZMIN,'O2_ZMIN',2,                    &
                          long_name='Vertical Minimum of O2',          &
                          units='mmol/m^3', grid_loc='2111',           &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_O2_ZMIN_DEPTH,'O2_ZMIN_DEPTH',2,        &
                          long_name='Depth of Vertical Minimum of O2', &
                          units='cm', grid_loc='2111',                 &
                          coordinates='TLONG TLAT time')

!-----------------------------------------------------------------------
!  nonstandard 3D fields
!-----------------------------------------------------------------------

   call define_tavg_field(tavg_O2_PRODUCTION,'O2_PRODUCTION',3,        &
                          long_name='O2 Production',                   &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_O2_CONSUMPTION,'O2_CONSUMPTION',3,      &
                          long_name='O2 Consumption',                  &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_AOU,'AOU',3,                            &
                          long_name='Apparent O2 Utilization ',        &
                          units='mmol/m^3', grid_loc='3111',           &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_PO4_RESTORE,'PO4_RESTORE',3,            &
                          long_name='PO4 Restoring',                   &
                          units='mmol/m^3', grid_loc='3111',           &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_NO3_RESTORE,'NO3_RESTORE',3,            &
                          long_name='NO3 Restoring',                   &
                          units='mmol/m^3', grid_loc='3111',           &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_SiO3_RESTORE,'SiO3_RESTORE',3,          &
                          long_name='SiO3 Restoring',                  &
                          units='mmol/m^3', grid_loc='3111',           &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_PAR_avg,'PAR_avg',3,                    &
                          long_name='PAR Average over Model Cell',     &
                          units='w/m^2', grid_loc='3114',              &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_POC_FLUX_IN,'POC_FLUX_IN',3,            &
                          long_name='POC Flux into Cell',              &
                          units='mmol/m^3 cm/s', grid_loc='3111',      &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_POC_PROD,'POC_PROD',3,                  &
                          long_name='POC Production',                  &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_POC_REMIN,'POC_REMIN',3,                &
                          long_name='POC Remineralization',            &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_POC_ACCUM,'POC_ACCUM',3,                &
                          long_name='POC Accumulation',                &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_CaCO3_FLUX_IN,'CaCO3_FLUX_IN',3,        &
                          long_name='CaCO3 flux into cell',            &
                          units='mmol/m^3 cm/s', grid_loc='3111',      &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_CaCO3_PROD,'CaCO3_PROD',3,              &
                          long_name='CaCO3 Production',                &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_CaCO3_REMIN,'CaCO3_REMIN',3,            &
                          long_name='CaCO3 Remineralization',          &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_SiO2_FLUX_IN,'SiO2_FLUX_IN',3,          &
                          long_name='SiO2 Flux into Cell',             &
                          units='mmol/m^3 cm/s', grid_loc='3111',      &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_SiO2_PROD,'SiO2_PROD',3,                &
                          long_name='SiO2 Production',                 &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_SiO2_REMIN,'SiO2_REMIN',3,              &
                          long_name='SiO2 Remineralization',           &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_dust_FLUX_IN,'dust_FLUX_IN',3,          &
                          long_name='Dust Flux into Cell',             &
                          units='ng/s/m^2', grid_loc='3111',           &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_dust_REMIN,'dust_REMIN',3,              &
                          long_name='Dust Remineralization',           &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_P_iron_FLUX_IN,'P_iron_FLUX_IN',3,      &
                          long_name='P_iron Flux into Cell',           &
                          units='mmol/m^3 cm/s', grid_loc='3111',      &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_P_iron_PROD,'P_iron_PROD',3,            &
                          long_name='P_iron Production',               &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_P_iron_REMIN,'P_iron_REMIN',3,          &
                          long_name='P_iron Remineralization',         &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_graze_sp,'graze_sp',3,                  &
                          long_name='Small Phyto Grazing',             &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_graze_diat,'graze_diat',3,              &
                          long_name='Diatom Grazing',                  &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_graze_diaz,'graze_diaz',3,              &
                          long_name='Diazotroph Grazing',              &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_graze_TOT,'graze_TOT',3,                &
                          long_name='Total Grazing',                   &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_sp_loss,'sp_loss',3,                    &
                          long_name='Small Phyto Loss',                &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_diat_loss,'diat_loss',3,                &
                          long_name='Diatom Loss',                     &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_diaz_loss,'diaz_loss',3,                &
                          long_name='Diaz Loss',                       &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_zoo_loss,'zoo_loss',3,                  &
                          long_name='Zooplankton Loss',                &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_sp_agg,'sp_agg',3,                      &
                          long_name='Small Phyto Aggregate',           &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_diat_agg,'diat_agg',3,                  &
                          long_name='Diatom Aggregate',                &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_photoC_sp,'photoC_sp',3,                &
                          long_name='Small Phyto C Fixation',          &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_CaCO3_form,'CaCO3_form',3,              &
                          long_name='CaCO3 Formation',                 &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_photoC_diat,'photoC_diat',3,            &
                          long_name='Diatom C Fixation',               &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_photoC_diaz,'photoC_diaz',3,            &
                          long_name='Diaz C Fixation',                 &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_photoC_TOT,'photoC_TOT',3,              &
                          long_name='Total C Fixation',                &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_photoC_sp_zint,'photoC_sp_zint',2,      &
                          long_name='Small Phyto C Fixation Vertical Integral',&
                          units='mmol/m^3 cm/s', grid_loc='2111',      &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_CaCO3_form_zint,'CaCO3_form_zint',2,    &
                          long_name='CaCO3 Formation Vertical Integral',&
                          units='mmol/m^3 cm/s', grid_loc='2111',      &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_photoC_diat_zint,'photoC_diat_zint',2,  &
                          long_name='Diatom C Fixation Vertical Integral',&
                          units='mmol/m^3 cm/s', grid_loc='2111',      &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_photoC_diaz_zint,'photoC_diaz_zint',2,  &
                          long_name='Diaz C Fixation Vertical Integral',&
                          units='mmol/m^3 cm/s', grid_loc='2111',      &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_photoC_TOT_zint,'photoC_TOT_zint',2,    &
                          long_name='Total C Fixation Vertical Integral',&
                          units='mmol/m^3 cm/s', grid_loc='2111',      &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_photoC_NO3_sp,'photoC_NO3_sp',3,        &
                          long_name='Small Phyto C Fixation from NO3',      &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_photoC_NO3_sp_zint,'photoC_NO3_sp_zint',2,&
                          long_name='Small Phyto C Fixation from NO3 Vertical Integral',&
                          units='mmol/m^3 cm/s', grid_loc='2111',      &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_photoC_NO3_diat,'photoC_NO3_diat',3,    &
                          long_name='Diatom C Fixation from NO3',           &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_photoC_NO3_diat_zint,'photoC_NO3_diat_zint',2,&
                          long_name='Diatom C Fixation from NO3 Vertical Integral',&
                          units='mmol/m^3 cm/s', grid_loc='2111',      &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_photoC_NO3_diaz,'photoC_NO3_diaz',3,    &
                          long_name='Diaz C Fixation from NO3',             &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_photoC_NO3_diaz_zint,'photoC_NO3_diaz_zint',2,&
                          long_name='Diaz C Fixation from NO3 Vertical Integral',&
                          units='mmol/m^3 cm/s', grid_loc='2111',      &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_photoC_NO3_TOT,'photoC_NO3_TOT',3,      &
                          long_name='Total C Fixation from NO3',            &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_photoC_NO3_TOT_zint,'photoC_NO3_TOT_zint',2,&
                          long_name='Total C Fixation from NO3 Vertical Integral',&
                          units='mmol/m^3 cm/s', grid_loc='2111',      &
                          coordinates='TLONG TLAT time')

!-----------------------------------------------------------------------
!  MORE nonstandard 3D fields
!-----------------------------------------------------------------------

   call define_tavg_field(tavg_DOC_prod,'DOC_prod',3,                  &
                          long_name='DOC Production',                  &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_DOC_remin,'DOC_remin',3,                &
                          long_name='DOC Remineralization',            &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_DON_prod,'DON_prod',3,                  &
                          long_name='DON Production',                  &
                          units='mmol/m^3/s', grid_loc='3114',          &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_DON_remin,'DON_remin',3,                &
                          long_name='DON Remineralization',            &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_DOFe_prod,'DOFe_prod',3,                &
                          long_name='DOFe Production',                 &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_DOFe_remin,'DOFe_remin',3,              &
                          long_name='DOFe Remineralization',           &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_DOP_prod,'DOP_prod',3,                  &
                          long_name='DOP Production',                  &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_DOP_remin,'DOP_remin',3,                &
                          long_name='DOP Remineralization',            &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_Fe_scavenge,'Fe_scavenge',3,            &
                          long_name='Iron Scavenging',                 &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_Fe_scavenge_rate,'Fe_scavenge_rate',3,  &
                          long_name='Iron Scavenging Rate',            &
                          units='1/y', grid_loc='3111',                &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_sp_N_lim,'sp_N_lim',3,                  &
                          long_name='Small Phyto N Limitation',        &
                          units='none', grid_loc='3114',               &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_sp_Fe_lim,'sp_Fe_lim',3,                &
                          long_name='Small Phyto Fe Limitation',       &
                          units='none', grid_loc='3114',               &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_sp_PO4_lim,'sp_PO4_lim',3,              &
                          long_name='Small Phyto PO4 Limitation',      &
                          units='none', grid_loc='3114',               &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_sp_light_lim,'sp_light_lim',3,          &
                          long_name='Small Phyto Light Limitation',    &
                          units='none', grid_loc='3114',               &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_diat_N_lim,'diat_N_lim',3,              &
                          long_name='Diatom N Limitation',             &
                          units='none', grid_loc='3114',               &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_diat_Fe_lim,'diat_Fe_lim',3,            &
                          long_name='Diatom Fe Limitation',            &
                          units='none', grid_loc='3114',               &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_diat_PO4_lim,'diat_PO4_lim',3,          &
                          long_name='Diatom PO4 Limitation',           &
                          units='none', grid_loc='3114',               &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_diat_SiO3_lim,'diat_SiO3_lim',3,        &
                          long_name='Diatom SiO3 Limitation',          &
                          units='none', grid_loc='3114',               &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_diat_light_lim,'diat_light_lim',3,      &
                          long_name='Diatom Light Limitation',         &
                          units='none', grid_loc='3114',               &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_diaz_Fe_lim,'diaz_Fe_lim',3,            &
                          long_name='Diaz Fe Limitation',              &
                          units='none', grid_loc='3114',               &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_diaz_P_lim,'diaz_P_lim',3,              &
                          long_name='Diaz PO4 Limitation',             &
                          units='none', grid_loc='3114',               &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_diaz_light_lim,'diaz_light_lim',3,      &
                          long_name='Diaz LIGHT Limitation',           &
                          units='none', grid_loc='3114',               &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_diaz_Nfix,'diaz_Nfix',3,                &
                          long_name='Diaz N Fixation',                 &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_bSi_form,'bSi_form',3,                  &
                          long_name='Diatom Si Uptake',                &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_NITRIF,'NITRIF',3,                      &
                          long_name='Nitrification',                   &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_DENITRIF,'DENITRIF',3,                  &
                          long_name='Denitrification',                 &
                          units='mmol/m^3/s', grid_loc='3111',         &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_photoFe_sp,'photoFe_sp',3,              &
                          long_name='Small Phyto Fe Uptake',           &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_photoFe_diat,'photoFe_diat',3,          &
                          long_name='Diatom Fe Uptake',                &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_photoFe_diaz,'photoFe_diaz',3,          &
                          long_name='Diaz Fe Uptake',                  &
                          units='mmol/m^3/s', grid_loc='3114',         &
                          coordinates='TLONG TLAT z_t_150m time')

   call define_tavg_field(tavg_CO3,'CO3',3,                            &
                          long_name='Carbonate Ion Concentration',     &
                          units='mmol/m^3', grid_loc='3111',           &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_HCO3,'HCO3',3,                          &
                          long_name='Bicarbonate Ion Concentration',   &
                          units='mmol/m^3', grid_loc='3111',           &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_H2CO3,'H2CO3',3,                        &
                          long_name='Carbonic Acid Concentration',     &
                          units='mmol/m^3', grid_loc='3111',           &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_pH_3D,'pH_3D',3,                        &
                          long_name='pH',                              &
                          units='none', grid_loc='3111',               &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_co3_sat_calc,'co3_sat_calc',3,          &
                          long_name='CO3 concentration at calcite saturation', &
                          units='mmol/m^3', grid_loc='3111',           &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_zsatcalc,'zsatcalc',2,                  &
                          long_name='Calcite Saturation Depth',        &
                          units='cm', grid_loc='2110',                 &
                          coordinates='TLONG TLAT time')

   call define_tavg_field(tavg_co3_sat_arag,'co3_sat_arag',3,          &
                          long_name='CO3 concentration at aragonite saturation', &
                          units='mmol/m^3', grid_loc='3111',           &
                          coordinates='TLONG TLAT z_t time')

   call define_tavg_field(tavg_zsatarag,'zsatarag',2,                  &
                          long_name='Aragonite Saturation Depth',      &
                          units='cm', grid_loc='2110',                 &
                          coordinates='TLONG TLAT time')

!-----------------------------------------------------------------------
!EOC

 end subroutine ecosys_init_tavg

!***********************************************************************
!BOP
! !IROUTINE: ecosys_set_interior
! !INTERFACE:

 subroutine ecosys_set_interior(k, TEMP_OLD, TEMP_CUR, SALT_OLD, SALT_CUR, &
    TRACER_MODULE_OLD, TRACER_MODULE_CUR, DTRACER_MODULE, this_block)

! !DESCRIPTION:
!  Compute time derivatives for ecosystem state variables
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   integer (int_kind), intent(in) :: &
      k                   ! vertical level index

   real (r8), dimension(nx_block,ny_block), intent(in) :: &
      TEMP_OLD,          &! old potential temperature (C)
      TEMP_CUR,          &! current potential temperature (C)
      SALT_OLD,          &! old salinity (msu)
      SALT_CUR            ! current salinity (msu)

   real (r8), dimension(nx_block,ny_block,km,ecosys_tracer_cnt), intent(in) :: &
      TRACER_MODULE_OLD, &! old tracer values
      TRACER_MODULE_CUR   ! current tracer values

   type (block), intent(in) :: &
      this_block          ! block info for the current block

! !OUTPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block,ecosys_tracer_cnt), intent(out) :: &
      DTRACER_MODULE      ! computed source/sink terms

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   character(*), parameter :: &
      subname = 'ecosys_mod:ecosys_set_interior'

   real (r8), parameter :: &
      epsC      = 1.00e-8, & ! small C concentration (mmol C/m^3)
      epsTinv   = 3.17e-8, & ! small inverse time scale (1/year) (1/sec)
      epsnondim = 1.00e-6    ! small non-dimensional number (non-dim)

   type(sinking_particle), save :: &
      POC,            & ! base units = nmol C
      P_CaCO3,        & ! base units = nmol CaCO3
      P_SiO2,         & ! base units = nmol SiO2
      dust,           & ! base units = g
      P_iron            ! base units = nmol Fe

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic), save :: &
      QA_dust_def,    & ! incoming deficit in the QA(dust) POC flux
      ZSATCALC,       & ! Calcite Saturation Depth
      ZSATARAG,       & ! Aragonite Saturation Depth
      CO3_CALC_ANOM_km1,&! CO3 concentration above calcite saturation at k-1
      CO3_ARAG_ANOM_km1 ! CO3 concentration above aragonite saturation at k-1

   real (r8), dimension(nx_block,ny_block) :: &
      TEMP,           & ! local copy of model TEMP
      SALT,           & ! local copy of model SALT
      DIC_loc,        & ! local copy of model DIC
      ALK_loc,        & ! local copy of model ALK
      PO4_loc,        & ! local copy of model PO4
      NO3_loc,        & ! local copy of model NO3
      SiO3_loc,       & ! local copy of model SiO3
      NH4_loc,        & ! local copy of model NH4
      Fe_loc,         & ! local copy of model Fe
      O2_loc,         & ! local copy of model O2
      DOC_loc,        & ! local copy of model DOC
      spC_loc,        & ! local copy of model spC
      spChl_loc,      & ! local copy of model spChl
      spCaCO3_loc,    & ! local copy of model spCaCO3
      diatC_loc,      & ! local copy of model diatC
      diatChl_loc,    & ! local copy of model diatChl
      zooC_loc,       & ! local copy of model zooC
      spFe_loc,       & ! local copy of model spFe
      diatSi_loc,     & ! local copy of model diatSi
      diatFe_loc,     & ! local copy of model diatFe
      diazC_loc,      & ! local copy of model diazC
      diazChl_loc,    & ! local copy of model diazChl
      diazFe_loc,     & ! local copy of model diazFe
      DON_loc,        & ! local copy of model DON
      DOFe_loc,       & ! local copy of model DOFe
      DOP_loc           ! local copy of model DOP

   real (r8), dimension(nx_block,ny_block) :: &
      WORK1,WORK2,WORK3 ! temporaries

   real (r8) :: &
      C_loss_thres,   & ! bio-C threshold at which losses go to zero (mmol C/m^3)
      f_loss_thres      ! fraction of grazing loss reduction at depth

   real (r8), dimension(nx_block,ny_block) :: &
      PAR_in,         & ! photosynthetically available radiation (W/m^2)
      KPARdz,         & ! PAR adsorption coefficient (non-dim)
      PAR_avg,        & ! average PAR over mixed layer depth (W/m^2)
      DOC_prod,       & ! production of DOC (mmol C/m^3/sec)
      DOC_remin,      & ! remineralization of DOC (mmol C/m^3/sec)
      NITRIF,         & ! nitrification (NH4 -> NO3) (mmol N/m^3/sec)
      DENITRIF,       & ! nitrification (NO3 -> N2) (mmol N/m^3/sec)
      RESTORE           ! restoring terms for nutrients (mmol ./m^3/sec)

   real (r8), dimension(nx_block,ny_block) :: &
      z_umax,         & ! max. zoo growth rate on sp at local T (1/sec)
      diat_umax,      & ! max. zoo growth rate on diatoms at local T (1/sec)
      z_mort,         & ! zoo respiration loss, (1/sec/((mmol C/m3))
      C_loss_diaz,    & ! bio-C threshold at which losses go to zero (mmol C/m^3)
      z_mort2,        & ! zoo quad mort rate, tohigherlevels (1/sec/((mmol C/m3))
      diaz_umax         ! max. zoo growth rate on diazotrophs at local T (1/sec)


   real (r8), dimension(nx_block,ny_block) :: &
      thetaC_sp,      & ! local Chl/C ratio in small phyto (mg Chl/mmol C)
      thetaC_diat,    & ! local Chl/C ratio in diatoms (mg Chl/mmol C)
      QCaCO3,         & ! small phyto CaCO3/C ratio (mmol CaCO3/mmol C)
      Tfunc,          & ! temp response function GD98 (non-dim)
      VNO3_sp,        & ! small phyto NO3 uptake rate (non-dim)
      VNH4_sp,        & ! small phyto NH4 uptake rate (non-dim)
      VNtot_sp,       & ! small phyto total N uptake rate (non-dim)
      VFeC_sp,        & ! ??? small phyto C-specific iron uptake (non-dim)
      VPO4_sp,        & ! ??? small phyto C-specific po4 uptake (non-dim)
      f_nut,          & ! nut limitation factor, modifies C fixation (non-dim)
      PCmax,          & ! max value of PCphoto at temperature TEMP (1/sec)
      PCphoto_sp,     & ! small C-specific rate of photosynth. (1/sec)
      photoC_sp,      & ! small phyto C-fixation (mmol C/m^3/sec)
      NO3_V_sp,       & ! nitrate uptake by small phyto (mmol NO3/m^3/sec)
      NH4_V_sp,       & ! ammonium uptake by small phyto (mmol NH4/m^3/sec)
      VNC_sp,         & ! small phyto C-specific N uptake rate (mmol N/mmol C/sec)
      pChl,           & ! Chl synth. regulation term (mg Chl/mmol N)
      photoacc_sp,    & ! Chl synth. term in photoadapt. (GD98) (mg Chl/m^3/sec)
      CaCO3_PROD,     & ! prod. of CaCO3 by small phyto (mmol CaCO3/m^3/sec)
      VNO3_diat,      & ! diatom nitrate uptake rate (non-dim)
      VNH4_diat,      & ! diatom ammonium uptake rate (non-dim)
      VNtot_diat,     & ! diatom total N uptake rate (non-dim)
      VFeC_diat,      & ! diatom C-specific iron uptake (non-dim)
      VPO4_diat,      & ! diatom C-specific PO4 uptake (non-dim)
      VSiO3_diat,     & ! C-specific SiO3 uptake for diatoms (non-dim)
      PCphoto_diat,   & ! diatom C-specific rate of photosynth. (1/sec)
      photoC_diat,    & ! diatom C-fixation (mmol C/m^3/sec)
      NO3_V_diat,     & ! nitrate uptake by diatoms (mmol NO3/m^3/sec)
      NH4_V_diat,     & ! ammonium uptake by diatoms (mmol NH4/m^3/sec)
      VNC_diat,       & ! diatom C-specific N uptake rate (mmol N/mmol C/sec)
      photoacc_diat,  & ! Chl synth. term in photoadapt. (GD98) (mg Chl/m^3/sec)
      reduceV,        & ! factor in nutrient uptake (mmol C/m^3)^2
      graze_sp,       & ! grazing rate on small phyto (mmol C/m^3/sec)
      graze_sp_zoo,   & ! graze_sp routed to zoo (mmol C/m^3/sec)
      graze_sp_poc,   & ! graze_sp routed to poc (mmol C/m^3/sec)
      graze_sp_doc,   & ! graze_sp routed to doc (mmol C/m^3/sec)
      graze_sp_dic      ! graze_sp routed to dic (mmol C/m^3/sec)

   real (r8), dimension(nx_block,ny_block) :: & ! max of 39 continuation lines
      graze_diat,     & ! grazing rate on diatoms (mmol C/m^3/sec)
      graze_diat_zoo, & ! graze_diat routed to zoo (mmol C/m^3/sec)
      graze_diat_poc, & ! graze_diat routed to poc (mmol C/m^3/sec)
      graze_diat_doc, & ! graze_diat routed to doc (mmol C/m^3/sec)
      graze_diat_dic, & ! graze_diat routed to dic (mmol C/m^3/sec)
      Pprime,         & ! used to limit phyto mort at low biomass (mmol C/m^3)
      sp_loss,        & ! small phyto non-grazing mort (mmol C/m^3/sec)
      sp_loss_poc,    & ! sp_loss routed to poc (mmol C/m^3/sec)
      sp_loss_doc,    & ! sp_loss routed to doc (mmol C/m^3/sec)
      sp_loss_dic,    & ! sp_loss routed to dic (mmol C/m^3/sec)
      sp_agg,         & ! small phyto agg loss (mmol C/m^3/sec)
      diat_loss,      & ! diatom non-grazing mort (mmol C/m^3/sec)
      diat_loss_poc,  & ! diat_loss routed to poc (mmol C/m^3/sec)
      diat_loss_doc,  & ! diat_loss routed to doc (mmol C/m^3/sec)
      diat_loss_dic,  & ! diat_loss routed to dic (mmol C/m^3/sec)
      diat_agg,       & ! diatom agg (mmol C/m^3/sec)
      f_zoo_detr,     & ! frac of zoo losses into large detrital pool (non-dim)
      Fe_scavenge,    & ! loss of dissolved iron, scavenging (mmol Fe/m^3/sec)
      Zprime,         & ! used to limit zoo mort at low biomass (mmol C/m^3)
      zoo_loss,       & ! mortality & higher trophic grazing on zooplankton (mmol C/m^3/sec)
      zoo_loss_doc,   & ! zoo_loss routed to doc (mmol C/m^3/sec)
      zoo_loss_dic,   & ! zoo_loss routed to dic (mmol C/m^3/sec)
      light_lim,      & ! light limitation factor
      Qsi,            & ! Diatom initial Si/C ratio (mmol Si/mmol C)
      gQsi,           & ! diatom Si/C ratio for growth (new biomass)
      Qfe_sp,         & ! small phyto init fe/C ratio (mmolFe/mmolC)
      gQfe_sp,        & ! small phyto fe/C for growth
      Qfe_diat,       & ! diatom init fe/C ratio
      gQfe_diat,      & ! diatom fe/C ratio for growth
      Qfe_diaz,       & ! diazotrophs init fe/C ratio
      gQfe_diaz         ! diazotroph fe/C ratio for new growth

   real (r8), dimension(nx_block,ny_block) :: & ! max of 39 continuation lines
      PCphoto_diaz,   & ! diazotroph C-specific rate of photosynth. (1/sec)
      photoC_diaz,    & ! diazotroph C-fixation (mmol C/m^3/sec)
      Vfec_diaz,      & ! diazotroph C-specific iron uptake (non-dim)
      Vpo4_diaz,      & ! diazotroph C-specific po4 uptake (non-dim)
      photoacc_diaz,  & ! Chl synth. term in photoadapt. (GD98) (mg Chl/m^3/sec)
      Vnc_diaz,       & ! diazotroph C-specific N uptake rate (mmol N/mmol C/sec)
      diaz_Nexcrete,  & ! diazotroph fixed N excretion
      diaz_Nfix,      & ! diazotroph total Nitrogen fixation (mmol N/m^3/sec)
      thetaC_diaz,    & ! local Chl/C ratio in diazotrophs (mg Chl/mmol C)
      photoFe_diaz,   & ! iron uptake by diazotrophs (mmolFe/m^3/sec)
      photoFe_diat,   & ! iron uptake by diatoms
      photoFe_sp,     & ! iron uptake by small phyto
      photoSi_diat,   & ! silicon uptake by diatoms (mmol Si/m^3/sec)
      remaining_diazP,& ! used in routing P from diazotroph losses
      diaz_loss,      & ! diazotroph non-grazing mort (mmol C/m^3/sec)
      diaz_loss_doc,  & ! mortality routed to DOM pool
      diaz_loss_poc,  & ! mortalitiy routed to POC pool
      diaz_loss_dic,  & ! mortality routed to remin
      diaz_loss_dop,  & ! P from mort routed to DOP pool
      diaz_loss_dip,  & ! P from mort routed to remin
                        ! P from diaz losses must be routed differently than
                        ! other elements to ensure that sinking detritus and
                        ! zooplankton pools get their fixed P/C ratios, the
                        ! remaining P is split evenly between DOP and PO4
      graze_diaz,     & ! grazing rate on diazotrophs (mmol C/m^3/sec)
      graze_diaz_poc, & ! grazing routed to sinking detr (mmol C/m^3/sec)
      graze_diaz_doc, & ! grazing routed to DOC (mmol C/m^3/sec)
      graze_diaz_dic, & ! grazing routed to remin (mmol C/m^3/sec)
      graze_diaz_zoo, & ! grazing routed to new zoo biomass
      DON_remin,      & ! portion of DON remineralized
      DOFe_remin,     & ! portion of DOFe remineralized
      DOP_remin,      & ! portion of DOP remineralized
      DOM_remin,      & ! fraction of DOM remineralized at current TEMP
      Fe_scavenge_rate  ! annual scavenging rate of iron as % of ambient

   real (r8), dimension(nx_block,ny_block) :: & ! max of 39 continuation lines
      DON_prod,       & ! production of dissolved organic N
      DOFe_prod,      & ! produciton of dissolved organic Fe
      DOP_prod,       & ! production of dissolved organic P
      po4_V_sp,       & ! sp uptake of po4 (mmol po4/m3/s)
      po4_V_diaz,     & ! diaz uptake of po4 (mmol po4/m3/s)
      dop_V_sp,       & ! sp uptake of dop (mmol dop/m3/s)
      dop_V_diaz,     & ! diaz uptake of dop (mmol dop/m3/s)
      cks,            & ! constant used in Fe quota modification
      cksi,           & ! constant used in Si quota modification
      VNO3_diaz,      & ! relative NO3 uptake by diazotrophs
      VNH4_diaz,      & ! relative NH4 uptake by diazotrophs
      VNtot_diaz,     & ! total relative N uptake by diazotrophs
      photoNO3_diaz,  & ! total nitrate uptake by diazotrophs
      photoNH4_diaz,  & ! total NH4 uptake by diazotrophs
      O2_PRODUCTION,  & ! O2 production
      O2_CONSUMPTION    ! O2 consumption

   real (r8), dimension(nx_block,ny_block) :: &
      CO3,            &! carbonate ion
      HCO3,           &! bicarbonate ion
      H2CO3,          &! carbonic acid
      OMEGA_CALC,     &! solubility ratio for aragonite
      OMEGA_ARAG       ! solubility ratio for calcite

   integer (int_kind) :: &
      bid,            & ! local_block id
      kk                ! index for looping over k levels

!-----------------------------------------------------------------------

   bid = this_block%local_id

   call timer_start(ecosys_interior_timer, block_id=bid)

   DTRACER_MODULE = c0

!-----------------------------------------------------------------------
!  exit immediately if computations are not to be performed
!-----------------------------------------------------------------------

   if (.not. lsource_sink) then
      call timer_stop(ecosys_interior_timer, block_id=bid)
      return
   endif

!-----------------------------------------------------------------------
!  create local copies of model tracers
!  treat negative values as zero
!  apply mask to local copies
!-----------------------------------------------------------------------

   TEMP         = p5*(TEMP_OLD + TEMP_CUR)
   SALT         = p5*(SALT_OLD + SALT_CUR)*salt_to_ppt

   DIC_loc      = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,dic_ind) + &
                              TRACER_MODULE_CUR(:,:,k,dic_ind)))
   ALK_loc      = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,alk_ind) + &
                              TRACER_MODULE_CUR(:,:,k,alk_ind)))
   PO4_loc      = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,po4_ind) + &
                              TRACER_MODULE_CUR(:,:,k,po4_ind)))
   NO3_loc      = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,no3_ind) + &
                              TRACER_MODULE_CUR(:,:,k,no3_ind)))
   SiO3_loc     = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,sio3_ind) + &
                              TRACER_MODULE_CUR(:,:,k,sio3_ind)))
   NH4_loc      = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,nh4_ind) + &
                              TRACER_MODULE_CUR(:,:,k,nh4_ind)))
   Fe_loc       = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,fe_ind) + &
                              TRACER_MODULE_CUR(:,:,k,fe_ind)))
   O2_loc       = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,o2_ind) + &
                              TRACER_MODULE_CUR(:,:,k,o2_ind)))
   DOC_loc      = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,doc_ind) + &
                              TRACER_MODULE_CUR(:,:,k,doc_ind)))
   spC_loc      = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,spC_ind) + &
                              TRACER_MODULE_CUR(:,:,k,spC_ind)))
   spChl_loc    = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,spChl_ind) + &
                              TRACER_MODULE_CUR(:,:,k,spChl_ind)))
   spCaCO3_loc  = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,spCaCO3_ind) + &
                              TRACER_MODULE_CUR(:,:,k,spCaCO3_ind)))
   diatC_loc    = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,diatC_ind) + &
                              TRACER_MODULE_CUR(:,:,k,diatC_ind)))
   diatChl_loc  = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,diatChl_ind) + &
                              TRACER_MODULE_CUR(:,:,k,diatChl_ind)))
   zooC_loc     = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,zooC_ind) + &
                              TRACER_MODULE_CUR(:,:,k,zooC_ind)))
   spFe_loc     = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,spFe_ind) + &
                              TRACER_MODULE_CUR(:,:,k,spFe_ind)))
   diatSi_loc   = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,diatSi_ind) + &
                              TRACER_MODULE_CUR(:,:,k,diatSi_ind)))
   diatFe_loc   = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,diatFe_ind) + &
                              TRACER_MODULE_CUR(:,:,k,diatFe_ind)))
   diazC_loc    = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,diazC_ind) + &
                              TRACER_MODULE_CUR(:,:,k,diazC_ind)))
   diazChl_loc  = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,diazChl_ind) + &
                              TRACER_MODULE_CUR(:,:,k,diazChl_ind)))
   diazFe_loc   = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,diazFe_ind) + &
                              TRACER_MODULE_CUR(:,:,k,diazFe_ind)))
   DON_loc      = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,don_ind) + &
                              TRACER_MODULE_CUR(:,:,k,don_ind)))
   DOFe_loc     = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,dofe_ind) + &
                              TRACER_MODULE_CUR(:,:,k,dofe_ind)))
   DOP_loc      = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,dop_ind) + &
                              TRACER_MODULE_CUR(:,:,k,dop_ind)))

   where (.not. LAND_MASK(:,:,bid) .or. k > KMT(:,:,bid))
      PO4_loc      = c0
      NO3_loc      = c0
      SiO3_loc     = c0
      NH4_loc      = c0
      Fe_loc       = c0
      O2_loc       = c0
      DOC_loc      = c0
      spC_loc      = c0
      spChl_loc    = c0
      spCaCO3_loc  = c0
      diatC_loc    = c0
      diatChl_loc  = c0
      zooC_loc     = c0
      spFe_loc     = c0
      diatSi_loc   = c0
      diatFe_loc   = c0
      diazC_loc    = c0
      diazChl_loc  = c0
      diazFe_loc   = c0
      DON_loc      = c0
      DOFe_loc     = c0
      DOP_loc      = c0
   end where

!-----------------------------------------------------------------------
!  If any phyto box are zero, set others to zeros.
!-----------------------------------------------------------------------

   where (spC_loc == c0 .or. spChl_loc == c0 .or. spFe_loc == c0)
      spC_loc     = c0
      spChl_loc   = c0
      spCaCO3_loc = c0
      spFe_loc    = c0
   end where

   where (diatC_loc == c0 .or. diatChl_loc == c0 .or. diatFe_loc == c0 &
             .or. diatSi_loc == c0)
      diatC_loc   = c0
      diatChl_loc = c0
      diatFe_loc  = c0
      diatSi_loc  = c0
   end where

   where (diazC_loc == c0 .or. diazChl_loc == c0 .or. diazFe_loc == c0)
      diazC_loc   = c0
      diazChl_loc = c0
      diazFe_loc  = c0
   end where

!-----------------------------------------------------------------------
!  set local variables, with incoming ratios
!-----------------------------------------------------------------------

   thetaC_sp   = spChl_loc   / (spC_loc + epsC)
   thetaC_diat = diatChl_loc / (diatC_loc + epsC)
   thetaC_diaz = diazChl_loc / (diazC_loc + epsC)
   Qsi         = diatSi_loc  / (diatC_loc + epsC)
   Qfe_diat    = diatFe_loc  / (diatC_loc + epsC)
   Qfe_sp      = spFe_loc    / (spC_loc + epsC)
   Qfe_diaz    = diazFe_loc  / (diazC_loc + epsC)
   where (Qsi > gQsi_max) Qsi = gQsi_max

!-----------------------------------------------------------------------
!  DETERMINE NEW ELEMENTAL RATIOS FOR GROWTH (NEW BIOMASS)
!-----------------------------------------------------------------------

   gQsi      = gQsi_0
   gQfe_diat = gQfe_diat_0
   gQfe_sp   = gQfe_sp_0
   gQfe_diaz = gQfe_diaz_0
   cks       = 10._r8
   cksi      = 10._r8

!-----------------------------------------------------------------------
!  Modify these initial ratios under low ambient iron conditions
!-----------------------------------------------------------------------

   where (Fe_loc < cks * parm_diat_kfe)
      gQfe_diat = max((gQfe_diat * Fe_loc /(cks * parm_diat_kfe)), &
                      gQfe_diat_min)
   end where

   where ((Fe_loc < cksi * parm_diat_kfe) .and. (Fe_loc > c0) .and. &
          (SiO3_loc > (cksi * parm_diat_kSiO3)))
      gQsi = min((gQsi*cksi*parm_diat_kfe/Fe_loc), gQsi_max)
   end where

   where (Fe_loc == c0)
      gQsi = gQsi_max
   end where

   where (Fe_loc < cks * parm_sp_kfe)
      gQfe_sp = max((gQfe_sp*Fe_loc/(cks * parm_sp_kfe)), &
                    gQfe_sp_min)
   end where

   where (Fe_loc < cks * diaz_kFe)
      gQfe_diaz = max((gQfe_diaz*Fe_loc/(cks * diaz_kFe)), &
                      gQfe_diaz_min)
   end where

!-----------------------------------------------------------------------
!  Modify the initial si/C ratio under low ambient Si conditions
!-----------------------------------------------------------------------

   where (SiO3_loc < (cksi * parm_diat_kSiO3))
      gQsi = max((gQsi*SiO3_loc/(cksi * parm_diat_kSiO3)), &
                 gQsi_min)
   end where

!-----------------------------------------------------------------------
!  various k==1 initializations
!
!  0.45   fraction of incoming SW -> PAR (non-dim)
!-----------------------------------------------------------------------

   if (k == 1) then
      call init_particulate_terms(POC, P_CaCO3, P_SiO2, dust, P_iron, &
                                  QA_dust_def, dust_FLUX_IN, this_block)
   endif

!-----------------------------------------------------------------------
!  QCaCO3 is the percentage of sp organic matter which is associated
!  with coccolithophores
!-----------------------------------------------------------------------

   QCaCO3 = spCaCO3_loc / (spC_loc + epsC)
   where (QCaCO3 > QCaCO3_max) QCaCO3 = QCaCO3_max

!-----------------------------------------------------------------------
!  compute PAR related quantities
!
!  0.03e-2   atten. coeff. per unit chlorophyll (1/cm/(mg Chl/m^3))
!  0.04e-2   atten. coeff. for water (1/cm)
!-----------------------------------------------------------------------

   PAR_in = PAR_out(:,:,bid)
   where (.not. LAND_MASK(:,:,bid) .or. k > KMT(:,:,bid)) PAR_in = c0

   if (partial_bottom_cells) then
      KPARdz = (k_chl * (spChl_loc + diatChl_loc + &
                         diazChl_loc) + k_h2o) * DZT(:,:,k,bid)
   else
      KPARdz = (k_chl * (spChl_loc + diatChl_loc + &
                         diazChl_loc) + k_h2o) * dz(k)
   endif

   PAR_out(:,:,bid) = PAR_in * exp(-KPARdz)
   PAR_avg = PAR_in * (c1 - exp(-KPARdz)) / KPARdz

!-----------------------------------------------------------------------
!  Tref = 30.0 reference temperature (deg. C)
!
!  Using q10 formulation with Q10 value of 2.0 (Doney et al., 1996).
!-----------------------------------------------------------------------

   Tfunc = Q_10**(((TEMP + T0_Kelvin)-(Tref + T0_Kelvin))/c10)

!-----------------------------------------------------------------------
!  modify growth mort rates by Tfunc
!-----------------------------------------------------------------------

   z_umax    = parm_z_umax_0    * Tfunc
   diat_umax = parm_diat_umax_0 * Tfunc
   z_mort2   = parm_z_mort2_0   * Tfunc
   z_mort    = parm_z_mort_0    * Tfunc
   diaz_umax = parm_diaz_umax_0 * Tfunc

   DOM_remin = parm_sd_remin_0

!-----------------------------------------------------------------------
!  Get relative nutrient uptake rates for phytoplankton,
!  min. relative uptake rate modifies C fixation in the manner
!  that the min. cell quota does in GD98.
!-----------------------------------------------------------------------

   VNO3_sp = (NO3_loc / parm_sp_kNO3) / &
      (c1 + (NO3_loc / parm_sp_kNO3) + (NH4_loc / parm_sp_kNH4))

   VNH4_sp = (NH4_loc / parm_sp_kNH4) / &
      (c1 + (NO3_loc / parm_sp_kNO3) + (NH4_loc / parm_sp_kNH4))

   VNtot_sp = VNO3_sp + VNH4_sp

   call accumulate_tavg_field(VNtot_sp, tavg_sp_N_lim,bid,k)

!-----------------------------------------------------------------------
!  get relative Fe uptake by phytoplankton
!  get relative P uptake rates for phytoplankton
!-----------------------------------------------------------------------

   VFeC_sp = Fe_loc / (Fe_loc + parm_sp_kFe)

   call accumulate_tavg_field(VFeC_sp, tavg_sp_Fe_lim,bid,k)

!-----------------------------------------------------------------------
!  Partition diaz P uptake in same manner to no3/nh4 partition
!-----------------------------------------------------------------------

   VPO4_sp = PO4_loc / (PO4_loc + parm_sp_kPO4)

   call accumulate_tavg_field(VPO4_sp, tavg_sp_PO4_lim,bid,k)

!-----------------------------------------------------------------------
!  Small Phytoplankton C-fixation - given light and relative uptake rates
!  determine small phyto nutrient limitation factor for carbon fixation
!-----------------------------------------------------------------------

   f_nut = min(VNtot_sp, VFeC_sp)
   f_nut = min(f_nut, VPO4_sp)

!-----------------------------------------------------------------------
!  get small phyto photosynth. rate, phyto C biomass change, photoadapt
!-----------------------------------------------------------------------

   PCmax = PCref * f_nut * Tfunc

   light_lim = (c1 - exp((-c1 * parm_alphaChlsp * thetaC_sp * PAR_avg) / &
                         (PCmax + epsTinv)))
   PCphoto_sp = PCmax * light_lim

   call accumulate_tavg_field(light_lim, tavg_sp_light_lim,bid,k)

   photoC_sp = PCphoto_sp * spC_loc

!-----------------------------------------------------------------------
!  Get nutrient uptakes by small phyto based on calculated C fixation
!  total N uptake (VNC_sp) is used in photoadaption
!-----------------------------------------------------------------------

   where (VNtot_sp > c0)
      NO3_V_sp = (VNO3_sp / VNtot_sp) * photoC_sp * Q
      NH4_V_sp = (VNH4_sp / VNtot_sp) * photoC_sp * Q
      VNC_sp = PCphoto_sp * Q
   elsewhere
      NO3_V_sp = c0
      NH4_V_sp = c0
      VNC_sp = c0
      photoC_sp = c0
   end where

   po4_V_sp = photoC_sp * Qp
   photoFe_sp = photoC_sp * gQfe_sp

   call accumulate_tavg_field(photoFe_sp, tavg_photoFe_sp,bid,k)

!-----------------------------------------------------------------------
!  calculate pChl, (used in photoadapt., GD98)
!  2.3   max value of thetaN (Chl/N ratio) (mg Chl/mmol N)
!  GD 98 Chl. synth. term
!-----------------------------------------------------------------------

   WORK1 = parm_alphaChlsp * thetaC_sp * PAR_avg
   where (WORK1 > c0)
      pChl = thetaN_max_sp * PCphoto_sp / WORK1
      photoacc_sp = (pChl * VNC_sp / thetaC_sp) * spChl_loc
   elsewhere
      photoacc_sp = c0
   end where

!-----------------------------------------------------------------------
!  CaCO3 Production, parameterized as function of small phyto production
!  decrease CaCO3 as function of nutrient limitation decrease CaCO3 prod
!  at low temperatures increase CaCO3 prod under bloom conditions
!  maximum calcification rate is 40% of primary production
!-----------------------------------------------------------------------

   CaCO3_PROD = f_prod_sp_CaCO3 * photoC_sp
   CaCO3_PROD = CaCO3_PROD * f_nut

   where (TEMP < CaCO3_temp_thres1)  &
      CaCO3_PROD = CaCO3_PROD * max((TEMP-CaCO3_temp_thres2), c0) / &
                   (CaCO3_temp_thres1-CaCO3_temp_thres2)

   where (spC_loc > CaCO3_sp_thres)  &
      CaCO3_PROD = min((CaCO3_PROD*spC_loc/CaCO3_sp_thres), &
                       (f_photosp_CaCO3*photoC_sp))

   call accumulate_tavg_field(CaCO3_PROD, tavg_CaCO3_form,bid,k)

   if (accumulate_tavg_now(tavg_CaCO3_form_zint)) then
      if (partial_bottom_cells) then
         WORK1 = DZT(:,:,k,bid) * CaCO3_PROD
      else
         WORK1 = dz(k) * CaCO3_PROD
      endif
      call accumulate_tavg_field(WORK1, tavg_CaCO3_form_zint,bid,k)
   endif

!-----------------------------------------------------------------------
!  Relative uptake rates for diatoms nitrate is VNO3, ammonium is VNH4
!-----------------------------------------------------------------------

   VNO3_diat = (NO3_loc / parm_diat_kNO3) / &
      (c1 + (NO3_loc / parm_diat_kNO3) + (NH4_loc / parm_diat_kNH4))

   VNH4_diat = (NH4_loc / parm_diat_kNH4) / &
      (c1 + (NO3_loc / parm_diat_kNO3) + (NH4_loc / parm_diat_kNH4))

   VNtot_diat = VNO3_diat + VNH4_diat

   call accumulate_tavg_field(VNtot_diat, tavg_diat_N_lim,bid,k)

!-----------------------------------------------------------------------
!  get relative Fe uptake by diatoms
!  get relative P uptake rates for diatoms
!  get relative SiO3 uptake rate for diatoms
!-----------------------------------------------------------------------

   VFeC_diat = Fe_loc / (Fe_loc + parm_diat_kFe)

   call accumulate_tavg_field(VFeC_diat, tavg_diat_Fe_lim,bid,k)

   VPO4_diat = PO4_loc / (PO4_loc + parm_diat_kPO4)

   call accumulate_tavg_field(VPO4_diat, tavg_diat_PO4_lim,bid,k)

   VSiO3_diat = SiO3_loc / (SiO3_loc + parm_diat_kSiO3)

   call accumulate_tavg_field(VSiO3_diat, tavg_diat_SiO3_lim,bid,k)

!-----------------------------------------------------------------------
!  Diatom carbon fixation and photoadapt.
!  determine diatom nutrient limitation factor for carbon fixation
!-----------------------------------------------------------------------

   f_nut = min(VNtot_diat, VFeC_diat)
   f_nut = min(f_nut, VSiO3_diat)
   f_nut = min(f_nut, VPO4_diat)

!-----------------------------------------------------------------------
!  get diatom photosynth. rate, phyto C biomass change, photoadapt
!-----------------------------------------------------------------------

   PCmax = PCref * f_nut * Tfunc

   light_lim = (c1 - exp((-c1 * parm_alphaChl * thetaC_diat * PAR_avg) / &
                         (PCmax + epsTinv)))
   PCphoto_diat = PCmax * light_lim

   call accumulate_tavg_field(light_lim, tavg_diat_light_lim,bid,k)

   photoC_diat = PCphoto_diat * diatC_loc

!-----------------------------------------------------------------------
!  Get nutrient uptake by diatoms based on C fixation
!-----------------------------------------------------------------------

   where (VNtot_diat > c0)
      NO3_V_diat = (VNO3_diat / VNtot_diat) * photoC_diat * Q
      NH4_V_diat = (VNH4_diat / VNtot_diat) * photoC_diat * Q
      VNC_diat = PCphoto_diat * Q
   elsewhere
      NO3_V_diat = c0
      NH4_V_diat = c0
      VNC_diat = c0
      photoC_diat = c0
   end where

   photoFe_diat = photoC_diat * gQfe_diat
   photoSi_diat = photoC_diat * gQsi

   call accumulate_tavg_field(photoFe_diat, tavg_photoFe_diat,bid,k)

   call accumulate_tavg_field(photoSi_diat, tavg_bSi_form,bid,k)

!-----------------------------------------------------------------------
!  calculate pChl, (used in photoadapt., GD98)
!  3.0   max value of thetaN (Chl/N ratio) (mg Chl/mmol N)
!  GD 98 Chl. synth. term
!-----------------------------------------------------------------------

   WORK1 = parm_alphaChl * thetaC_diat * PAR_avg
   where (WORK1 > c0)
      pChl = thetaN_max_diat * PCphoto_diat / WORK1
      photoacc_diat = (pChl * VNC_diat / thetaC_diat) * diatChl_loc
   elsewhere
      photoacc_diat = c0
   end where

!-----------------------------------------------------------------------
!  get relative Fe uptake by diazotrophs
!  get relative P uptake rates for diazotrophs
!-----------------------------------------------------------------------

   Vfec_diaz = Fe_loc/(Fe_loc + diaz_kFe)

   call accumulate_tavg_field(Vfec_diaz, tavg_diaz_Fe_lim,bid,k)

   Vpo4_diaz = PO4_loc / (PO4_loc + diaz_kPO4)

   call accumulate_tavg_field(Vpo4_diaz, tavg_diaz_P_lim,bid,k)

   f_nut = min(Vpo4_diaz, Vfec_diaz)

!-----------------------------------------------------------------------
!  get diazotroph photosynth. rate, phyto C biomass change
!-----------------------------------------------------------------------

   PCmax = PCrefDiaz * f_nut * Tfunc
   where (TEMP < diaz_temp_thres) PCmax = c0

   light_lim = (c1 - exp((-c1 * parm_alphaDiaz * thetaC_diaz * PAR_avg) / &
                         (PCmax + epsTinv)))
   PCphoto_diaz = PCmax * light_lim

   call accumulate_tavg_field(light_lim, tavg_diaz_light_lim,bid,k)

   photoC_diaz = PCphoto_diaz * diazC_loc

!-----------------------------------------------------------------------
!  Get Fe and po4 uptake by diazotrophs based on C fixation
!-----------------------------------------------------------------------

   po4_V_diaz = photoC_diaz * Qp_diaz
   photoFe_diaz = photoC_diaz * gQfe_diaz

   call accumulate_tavg_field(photoFe_diaz, tavg_photoFe_diaz,bid,k)

!-----------------------------------------------------------------------
!  Relative uptake rates for diazotrophs nitrate is VNO3, ammonium is VNH4
!-----------------------------------------------------------------------

   VNO3_diaz = (NO3_loc / diaz_kNO3) / &
      (c1 + (NO3_loc / diaz_kNO3) + (NH4_loc / diaz_kNH4))

   VNH4_diaz = (NH4_loc / diaz_kNH4) / &
      (c1 + (NO3_loc / diaz_kNO3) + (NH4_loc / diaz_kNH4))

   VNtot_diaz = c1

!------------------------------------------------------------------------
!  Compute inoranic N uptake by diazotrophs.
!------------------------------------------------------------------------

   WORK1 = photoC_diaz * Q
   photoNO3_diaz = VNO3_diaz / VNtot_diaz *WORK1
   photoNH4_diaz = VNH4_diaz / VNtot_diaz *WORK1

!-----------------------------------------------------------------------
!  Get N fixation by diazotrophs based on C fixation,
!  Diazotrophs fix more than they need then 30% is excreted
!-----------------------------------------------------------------------

   diaz_Nfix     = (WORK1 * r_Nfix_photo) - photoNO3_diaz - photoNH4_diaz
   diaz_Nexcrete = diaz_Nfix + photoNO3_diaz + photoNH4_diaz - WORK1

   Vnc_diaz = PCphoto_diaz * Q

   call accumulate_tavg_field(diaz_Nfix, tavg_diaz_Nfix,bid,k)

!-----------------------------------------------------------------------
!  calculate pChl, (used in photoadapt., GD98)
!  3.4   max value of thetaN (Chl/N ratio) (mg Chl/mmol N)
!  GD 98 Chl. synth. term
!-----------------------------------------------------------------------

   WORK1 = parm_alphaDiaz * thetaC_diaz * PAR_avg
   where (WORK1 > c0)
      pChl = thetaN_max_diaz * PCphoto_diaz / WORK1
      photoacc_diaz = (pChl * Vnc_diaz / thetaC_diaz) * diazChl_loc
   elsewhere
      photoacc_diaz = c0
   end where

!-----------------------------------------------------------------------
!  CALCULATE GRAZING AND OTHER MORT
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!  calculate the loss threshold interpolation factor
!-----------------------------------------------------------------------

   if (zt(k) > thres_z1) then
      if (zt(k) < thres_z2) then
         f_loss_thres = (thres_z2 - zt(k))/(thres_z2 - thres_z1)
      else
         f_loss_thres = c0
      endif
   else
      f_loss_thres = c1
   endif

!-----------------------------------------------------------------------
!  0.001 small phytoplankton threshold C concentration (mmol C/m^3)
!  get small phyto loss(in C units)
!  small phyto agg loss
!  get grazing rate (graze_sp) on small phyto  (in C units)
!-----------------------------------------------------------------------

   C_loss_thres = f_loss_thres * loss_thres_sp

   Pprime = max(spC_loc - C_loss_thres, c0)

   sp_loss = sp_mort * Pprime * Tfunc

   sp_agg = min((sp_agg_rate_max * dps) * Pprime, sp_mort2 * Pprime * Pprime)

   reduceV = Pprime * Pprime
   graze_sp = z_umax * zooC_loc * (reduceV / (reduceV + parm_z_grz))

!-----------------------------------------------------------------------
!  routing of graze_sp & sp_loss
!  sp_agg all goes to POC
!  currently assumes that 33% of grazed caco3 is remineralized
!  if z_ingest ever changes, coefficients on routing grazed sp must change!
!  min.%C routed to POC from grazing for ballast requirements = 0.4 * Qcaco3
!  min.%C routed from sp_loss = 0.59 * QCaCO3, or P_CaCO3%rho
!-----------------------------------------------------------------------

   graze_sp_zoo = z_ingest * graze_sp
   graze_sp_poc = graze_sp * max((caco3_poc_min * QCaCO3),  &
                                 min(max(0.05_r8,(spc_poc_fac * Pprime)),&
                                     f_graze_sp_poc_lim))

   graze_sp_doc = f_graze_sp_doc * graze_sp
   graze_sp_dic = f_graze_sp_dic * graze_sp - graze_sp_poc

   sp_loss_poc = QCaCO3 * sp_loss
   sp_loss_doc = (c1 - parm_labile_ratio) * (sp_loss - sp_loss_poc)
   sp_loss_dic = parm_labile_ratio * (sp_loss - sp_loss_poc)

!-----------------------------------------------------------------------
!  0.01 small diatom threshold C concentration (mmol C/m^3)
!  get diatom loss(in C units)
!  Diatom agg loss, min. 1%/day
!  get grazing rate (graze_diat) on diatoms  (in C units)
!-----------------------------------------------------------------------

   C_loss_thres = f_loss_thres * loss_thres_diat

   Pprime = max(diatC_loc - C_loss_thres, c0)

   diat_loss = diat_mort * Pprime * Tfunc

   diat_agg = min((diat_agg_rate_max * dps) * Pprime, diat_mort2 * Pprime * Pprime)
   diat_agg = max((diat_agg_rate_min * dps) * Pprime, diat_agg)

!-----------------------------------------------------------------------
!  Lower z_grz term for diatoms and diazotrophs, larger, more mobile predators
!-----------------------------------------------------------------------

   reduceV = Pprime * Pprime
   graze_diat = diat_umax * zooC_loc * &
      (reduceV / (reduceV + 0.7_r8))

!-----------------------------------------------------------------------
!  routing of graze_diat & diat_loss
!  diat_agg all goes to POC
!  NOTE: if z_ingest is changed, coeff.s for poc,doc and dic must change!
!-----------------------------------------------------------------------

   graze_diat_zoo = z_ingest * graze_diat
   graze_diat_poc = f_graze_diat_poc * graze_diat
   graze_diat_doc = f_graze_diat_doc * graze_diat
   graze_diat_dic = f_graze_diat_dic * graze_diat

   diat_loss_poc = f_diat_loss_poc * diat_loss
   diat_loss_doc = (c1-parm_labile_ratio) * f_diat_loss_dc * diat_loss
   diat_loss_dic = parm_labile_ratio * f_diat_loss_dc * diat_loss

!-----------------------------------------------------------------------
!  0.03 small diazotroph threshold C concentration (mmol C/m^3)
!  Lower value used at temperatures < 16 deg. C, negligible biomass
!  get diazotroph loss(in C units)
!  get grazing rate (graze_diaz) on diazotrophs  (in C units)
!  no aggregation loss for diazotrophs
!-----------------------------------------------------------------------

   C_loss_diaz = f_loss_thres * loss_thres_diaz

   where (TEMP .lt. diaz_temp_thres) &
      C_loss_diaz = f_loss_thres * loss_thres_diaz2

   Pprime = max(diazC_loc - C_loss_diaz, c0)

   diaz_loss = diaz_mort * Pprime * Tfunc

   reduceV = Pprime * Pprime
   graze_diaz = diaz_umax * zooC_loc * (reduceV / (reduceV + parm_z_grz))

!-----------------------------------------------------------------------
!  routing of graze_diaz & diaz_loss
!-----------------------------------------------------------------------

   graze_diaz_zoo = f_graze_diaz_zoo * graze_diaz
   graze_diaz_poc = f_graze_diaz_poc * graze_diaz
   graze_diaz_doc = f_graze_diaz_doc * graze_diaz
   graze_diaz_dic = f_graze_diaz_dic * graze_diaz

   diaz_loss_poc = f_diaz_loss_poc * diaz_loss
   diaz_loss_doc = (c1 - parm_labile_ratio) * (diaz_loss - diaz_loss_poc)
   diaz_loss_dic = parm_labile_ratio * (diaz_loss - diaz_loss_poc)

!-----------------------------------------------------------------------
!  Note as diazotrophs have different Qp, we must route enough P into zoopl
!  and sinking detritus pools to fill their fixed p/C ratios.  The left over
!  P (remaining_diazP) is split between DOP and DIP pools
!-----------------------------------------------------------------------

   remaining_diazP =  ((graze_diaz + diaz_loss) * Qp_diaz) - &
                      ((graze_diaz_poc+graze_diaz_zoo) * Qp)
   diaz_loss_dop = (c1 - parm_labile_ratio) * remaining_diazP
   diaz_loss_dip = parm_labile_ratio * remaining_diazP

!-----------------------------------------------------------------------
!  get fractional factor for routing of zoo losses, based on food supply
!  more material is routed to large detrital pool when diatoms eaten
!-----------------------------------------------------------------------

   f_zoo_detr = (f_diat_zoo_detr * (graze_diat + epsC * epsTinv) + &
                 f_sp_zoo_detr * (graze_sp + epsC * epsTinv) + &
                 f_diaz_zoo_detr * (graze_diaz + epsC * epsTinv)) / &
                (graze_diat + graze_sp + graze_diaz + c3 * epsC * epsTinv)

!-----------------------------------------------------------------------
!  0.01 small zoo threshold C concentration (mmol C/m^3)
!  zoo losses
!-----------------------------------------------------------------------

   C_loss_thres = f_loss_thres * loss_thres_zoo

   Zprime = max(zooC_loc - C_loss_thres, c0)

   zoo_loss = z_mort2 * Zprime**1.4_r8 + z_mort * Zprime

   zoo_loss_doc = (c1 - parm_labile_ratio) * (c1 - f_zoo_detr) * zoo_loss
   zoo_loss_dic = parm_labile_ratio * (c1 - f_zoo_detr) * zoo_loss

!-----------------------------------------------------------------------
!  compute terms for DOM
!-----------------------------------------------------------------------

   DOC_prod = sp_loss_doc + graze_sp_doc + zoo_loss_doc + diat_loss_doc &
              + graze_diat_doc + diaz_loss_doc + graze_diaz_doc

   DON_prod = DOC_prod * Q
   DOP_prod = (sp_loss_doc + graze_sp_doc + zoo_loss_doc + diat_loss_doc &
               + graze_diat_doc) * Qp + diaz_loss_dop
   DOFe_prod = (zoo_loss_doc * Qfe_zoo) &
               + (Qfe_sp * (graze_sp_doc + sp_loss_doc)) &
               + (Qfe_diat * (graze_diat_doc + diat_loss_doc)) &
               + (Qfe_diaz * (graze_diaz_doc + diaz_loss_doc))

   DOC_remin = DOC_loc * DOM_remin
   DON_remin = DON_loc * DOM_remin
   DOFe_remin = DOFe_loc * DOM_remin
   DOP_remin = DOP_loc * DOM_remin

!-----------------------------------------------------------------------
!  large detritus C
!-----------------------------------------------------------------------

   POC%prod(:,:,bid) = sp_agg + graze_sp_poc + sp_loss_poc + &
                       f_zoo_detr * zoo_loss + diat_loss_poc + &
                       diat_agg + graze_diat_poc + graze_diaz_poc

!-----------------------------------------------------------------------
!  large detrital CaCO3
!  33% of CaCO3 is remin when phyto are grazed
!-----------------------------------------------------------------------

   P_CaCO3%prod(:,:,bid) = ((c1 - f_graze_CaCO3_REMIN) * graze_sp + &
                            sp_loss + sp_agg) * QCaCO3

!-----------------------------------------------------------------------
!  large detritus SiO2
!  grazed diatom SiO2, 60% is remineralized
!-----------------------------------------------------------------------

   P_SiO2%prod(:,:,bid) = ((c1 - f_graze_si_remin) * graze_diat + &
                           diat_agg  + f_diat_loss_poc * diat_loss) * Qsi

   dust%prod(:,:,bid) = c0

!-----------------------------------------------------------------------
!  Compute iron scavenging :
!  1) compute in terms of loss per year per unit iron (%/year/fe)
!  2) scale by sinking POMx6 + Dust + bSi + CaCO3 flux
!  3) increase scavenging at higher iron (>0.6nM)
!  4) convert to net loss per second
!-----------------------------------------------------------------------

   Fe_scavenge_rate = Fe_scavenge_rate0

   Fe_scavenge_rate = Fe_scavenge_rate * &
      ((POC%sflux_out(:,:,bid) + POC%hflux_out(:,:,bid)) * 72.06_r8 + &
       (P_CaCO3%sflux_out(:,:,bid) + P_CaCO3%hflux_out(:,:,bid)) * P_CaCO3%mass + &
       (P_SiO2%sflux_out(:,:,bid)+P_SiO2%hflux_out(:,:,bid))*P_SiO2%mass + &
       (dust%sflux_out(:,:,bid) + dust%hflux_out(:,:,bid)) * dust_fescav_scale)

   where (Fe_loc > Fe_scavenge_thres1) &
      Fe_scavenge_rate = Fe_scavenge_rate + &
                         (Fe_loc - Fe_scavenge_thres1) * fe_max_scale2

   Fe_scavenge = yps * Fe_loc * Fe_scavenge_rate

   P_iron%prod(:,:,bid) = &
      ((sp_agg + graze_sp_poc + sp_loss_poc) * Qfe_sp) &
      + (zoo_loss * f_zoo_detr * Qfe_zoo) &
      + ((diat_agg + graze_diat_poc + diat_loss_poc) * Qfe_diat) &
      + (graze_diaz_poc * Qfe_diaz) + (f_fescav_P_iron * Fe_scavenge)

   call compute_particulate_terms(k, POC, P_CaCO3, P_SiO2, dust, P_iron, &
                                  QA_dust_def, TEMP, O2_loc, this_block)

!-----------------------------------------------------------------------
!  nitrate & ammonium
!  nitrification in low light
!  use exponential decay of PAR across model level to compute taper factor
!-----------------------------------------------------------------------

   if (lrest_no3) then
      if (lnutr_variable_restore) then
         RESTORE = NUTR_RESTORE_RTAU(:,:,bid) * &
                      merge((NO3_CLIM(:,:,k,bid) - NO3_loc), &
                            c0, k <= NUTR_RESTORE_MAX_LEVEL(:,:,bid))
      else
         RESTORE = (NO3_CLIM(:,:,k,bid) - NO3_loc) * nutr_rest_time_inv(k)
      endif
   else
      RESTORE = c0
   endif

   call accumulate_tavg_field(RESTORE, tavg_NO3_RESTORE,bid,k)

   where (PAR_out(:,:,bid) < parm_nitrif_par_lim)
      NITRIF = parm_kappa_nitrif * NH4_loc
      where (PAR_in > parm_nitrif_par_lim)
         NITRIF = NITRIF * log(PAR_out(:,:,bid) / parm_nitrif_par_lim) / (-KPARdz)
      end where
   elsewhere
      NITRIF = c0
   end where

   call accumulate_tavg_field(NITRIF, tavg_NITRIF,bid,k)

!-----------------------------------------------------------------------
!  Compute denitrification under low O2 conditions
!-----------------------------------------------------------------------

   WORK1 = ((parm_o2_min+parm_o2_min_delta) - O2_loc) / parm_o2_min_delta
   WORK1 = min(max(WORK1,c0),c1)
   where (NO3_loc .lt. parm_no3_min)
      WORK1 = WORK1 * (NO3_loc / parm_no3_min)
   end where
   DENITRIF = WORK1 * (DOC_remin + POC%remin(:,:,bid)) / denitrif_C_N

   call accumulate_tavg_field(DENITRIF, tavg_DENITRIF,bid,k)

!-----------------------------------------------------------------------
!  nitrate & ammonium
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,no3_ind) = RESTORE + NITRIF - NO3_V_diat - &
      NO3_V_sp - DENITRIF - photoNO3_diaz

   DTRACER_MODULE(:,:,nh4_ind) = -NH4_V_diat - NH4_V_sp - NITRIF + &
        Q * (zoo_loss_dic + sp_loss_dic + graze_sp_dic + diat_loss_dic + &
        graze_diat_dic + POC%remin(:,:,bid) + diaz_loss_dic + &
        graze_diaz_dic) + DON_remin + diaz_Nexcrete - photoNH4_diaz

!-----------------------------------------------------------------------
!  dissolved iron
!-----------------------------------------------------------------------

     DTRACER_MODULE(:,:,fe_ind) = P_iron%remin(:,:,bid)  &
       + (Qfe_zoo * zoo_loss_dic) + DOFe_remin  &
       + (Qfe_sp * (sp_loss_dic + graze_sp_dic))     &
       + (Qfe_diat * (diat_loss_dic + graze_diat_dic))     &
       + (Qfe_diaz * (diaz_loss_dic + graze_diaz_dic))     &
       + graze_diaz_zoo * (Qfe_diaz-Qfe_zoo)     &
       + graze_diat_zoo * (Qfe_diat-Qfe_zoo)     &
       + graze_sp_zoo * (Qfe_sp-Qfe_zoo)     &
       - photoFe_sp - photoFe_diat - photoFe_diaz - Fe_scavenge

!-----------------------------------------------------------------------
!  dissolved SiO3
!-----------------------------------------------------------------------

   if (lrest_sio3) then
      if (lnutr_variable_restore) then
         RESTORE = NUTR_RESTORE_RTAU(:,:,bid) * &
                      merge((SiO3_CLIM(:,:,k,bid) - SiO3_loc), &
                            c0, k <= NUTR_RESTORE_MAX_LEVEL(:,:,bid))
      else
         RESTORE = (SiO3_CLIM(:,:,k,bid) - SiO3_loc) * nutr_rest_time_inv(k)
      endif
   else
      RESTORE = c0
   endif

   call accumulate_tavg_field(RESTORE, tavg_SiO3_RESTORE,bid,k)

   DTRACER_MODULE(:,:,sio3_ind) = RESTORE + P_SiO2%remin(:,:,bid) + &
      Qsi * (f_graze_si_remin * graze_diat +f_diat_loss_dc * diat_loss) &
      - photoSi_diat

!-----------------------------------------------------------------------
!  phosphate
!-----------------------------------------------------------------------

   if (lrest_po4) then
      if (lnutr_variable_restore) then
         RESTORE = NUTR_RESTORE_RTAU(:,:,bid) * &
                      merge((PO4_CLIM(:,:,k,bid) - PO4_loc), &
                            c0, k <= NUTR_RESTORE_MAX_LEVEL(:,:,bid))
      else
         RESTORE = (PO4_CLIM(:,:,k,bid) - PO4_loc) * nutr_rest_time_inv(k)
      endif
   else
      RESTORE = c0
   endif

   call accumulate_tavg_field(RESTORE, tavg_PO4_RESTORE,bid,k)

   DTRACER_MODULE(:,:,po4_ind) = RESTORE + (Qp * (POC%remin(:,:,bid) + &
      zoo_loss_dic + sp_loss_dic + graze_sp_dic + diat_loss_dic + &
      graze_diat_dic - photoC_diat)) &
      + DOP_remin + diaz_loss_dip - po4_V_sp - po4_V_diaz

!-----------------------------------------------------------------------
!  small phyto Carbon
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,spC_ind) = photoC_sp - graze_sp - sp_loss - sp_agg

!-----------------------------------------------------------------------
!  small phyto Chlorophyll
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,spChl_ind) = photoacc_sp - &
      thetaC_sp * (graze_sp + sp_loss + sp_agg)

!-----------------------------------------------------------------------
!  small phytoplankton CaCO3
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,spCaCO3_ind) = CaCO3_PROD - &
      (graze_sp + sp_loss + sp_agg) * QCaCO3

!-----------------------------------------------------------------------
!  diatom Carbon
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,diatC_ind) = photoC_diat - graze_diat - &
      diat_loss - diat_agg

!-----------------------------------------------------------------------
!  diatom Chlorophyll
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,diatChl_ind) = photoacc_diat - &
      thetaC_diat * (graze_diat + diat_loss + diat_agg)

!-----------------------------------------------------------------------
!  zoo Carbon
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,zooC_ind) = graze_sp_zoo + graze_diat_zoo &
      + graze_diaz_zoo - zoo_loss

!-----------------------------------------------------------------------
!  dissolved organic Matter
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,doc_ind) = DOC_prod - DOC_remin

   DTRACER_MODULE(:,:,don_ind) = DON_prod - DON_remin

   DTRACER_MODULE(:,:,dop_ind) = DOP_prod - DOP_remin

   DTRACER_MODULE(:,:,dofe_ind) = DOFe_prod - DOFe_remin

!-----------------------------------------------------------------------
!  small phyto Fe
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,spFe_ind) =  photoFe_sp &
      - (Qfe_sp * (graze_sp+sp_loss+sp_agg))

!-----------------------------------------------------------------------
!  Diatom Fe
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,diatFe_ind) =  photoFe_diat &
      - (Qfe_diat * (graze_diat+diat_loss+diat_agg))

!-----------------------------------------------------------------------
!  Diatom Si
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,diatSi_ind) =  photoSi_diat &
      - (Qsi * (graze_diat+diat_loss+diat_agg))

!-----------------------------------------------------------------------
!  Diazotroph C
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,diazC_ind) =  photoC_diaz - graze_diaz - diaz_loss

!-----------------------------------------------------------------------
!  diazotroph Chlorophyll
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,diazChl_ind) = photoacc_diaz - &
      thetaC_diaz * (graze_diaz + diaz_loss)

!-----------------------------------------------------------------------
!  Diazotroph Fe
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,diazFe_ind) =  photoFe_diaz &
      - (Qfe_diaz * (graze_diaz + diaz_loss))

!-----------------------------------------------------------------------
!   dissolved inorganic Carbon
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,dic_ind) = &
      DOC_remin + POC%remin(:,:,bid) + P_CaCO3%remin(:,:,bid) + &
      f_graze_CaCO3_REMIN * graze_sp * QCaCO3 + zoo_loss_dic + sp_loss_dic + &
      graze_sp_dic + diat_loss_dic + graze_diat_dic - photoC_sp - photoC_diat &
      - CaCO3_PROD + graze_diaz_dic + diaz_loss_dic - photoC_diaz

!-----------------------------------------------------------------------
!  alkalinity
!-----------------------------------------------------------------------

   DTRACER_MODULE(:,:,alk_ind) = -DTRACER_MODULE(:,:,no3_ind) + &
      DTRACER_MODULE(:,:,nh4_ind) + &
      c2 * (P_CaCO3%remin(:,:,bid) + f_graze_CaCO3_REMIN * graze_sp * QCaCO3 &
            - CaCO3_PROD)

!-----------------------------------------------------------------------
!  oxygen
!-----------------------------------------------------------------------

!  DTRACER_MODULE(:,:,o2_ind) = (photoC_sp + photoC_diat) / &
!     parm_Red_D_C_O2 + photoC_diaz/parm_Red_D_C_O2_diaz

!  WORK1 = (O2_loc - parm_o2_min) / parm_o2_min_delta
!  WORK1 = min(max(WORK1,c0),c1)
!  DTRACER_MODULE(:,:,o2_ind) = DTRACER_MODULE(:,:,o2_ind) + WORK1 &
!     * ((- POC%remin(:,:,bid) - DOC_remin - zoo_loss_dic - sp_loss_dic  &
!        - graze_sp_dic - diat_loss_dic - graze_diat_dic   &
!        - graze_diaz_dic - diaz_loss_dic) / parm_Red_P_C_O2)

   DTRACER_MODULE(:,:,o2_ind) = c0

   O2_PRODUCTION = c0

   where (photoC_diat > c0)
      O2_PRODUCTION = O2_PRODUCTION + photoC_diat * &
         ((NO3_V_diat/(NO3_V_diat+NH4_V_diat)) / parm_Red_D_C_O2 + &
          (NH4_V_diat/(NO3_V_diat+NH4_V_diat)) / parm_Remin_D_C_O2)
   end where

   where (photoC_sp > c0)
      O2_PRODUCTION = O2_PRODUCTION + photoC_sp * &
         ((NO3_V_sp/(NO3_V_sp+NH4_V_sp)) / parm_Red_D_C_O2 + &
          (NH4_V_sp/(NO3_V_sp+NH4_V_sp)) / parm_Remin_D_C_O2)
   end where

   where (photoC_diaz > c0)
      O2_PRODUCTION = O2_PRODUCTION + photoC_diaz * &
         ((photoNO3_diaz/(photoNO3_diaz+photoNH4_diaz+diaz_Nfix)) / parm_Red_D_C_O2 + &
          (photoNH4_diaz/(photoNO3_diaz+photoNH4_diaz+diaz_Nfix)) / parm_Remin_D_C_O2 + &
          (diaz_Nfix/(photoNO3_diaz+photoNH4_diaz+diaz_Nfix)) / parm_Red_D_C_O2_diaz)
   end where

   WORK1 = (O2_loc - parm_o2_min) / parm_o2_min_delta
   WORK1 = min(max(WORK1,c0),c1)
   O2_CONSUMPTION = WORK1 * &
      ((POC%remin(:,:,bid) + DOC_remin + zoo_loss_dic + &
        sp_loss_dic + graze_sp_dic + diat_loss_dic + graze_diat_dic + &
        graze_diaz_dic + diaz_loss_dic)/ parm_Remin_D_C_O2 + (c2*NITRIF))

   DTRACER_MODULE(:,:,o2_ind) = O2_PRODUCTION - O2_CONSUMPTION

!-----------------------------------------------------------------------
!  various tavg/history variables
!-----------------------------------------------------------------------

   if (k == 1) then
      if (accumulate_tavg_now(tavg_O2_ZMIN) .or. accumulate_tavg_now(tavg_O2_ZMIN_DEPTH)) then
         ! WORK1 = O2 at this level
         ! WORK2 = vertical min of O2
         ! WORK3 = depth of min

         kk = 1
         WORK1 = p5*(TRACER_MODULE_OLD(:,:,kk,o2_ind) + &
                     TRACER_MODULE_CUR(:,:,kk,o2_ind))
         WORK2 = WORK1
         WORK3 = zt(kk)

         do kk = 2,km
            WORK1 = p5*(TRACER_MODULE_OLD(:,:,kk,o2_ind) + &
                        TRACER_MODULE_CUR(:,:,kk,o2_ind))
            where (kk <= KMT(:,:,bid) .and. (WORK1 < WORK2))
               WORK2 = WORK1
               WORK3 = zt(kk)
            endwhere
         end do

         call accumulate_tavg_field(WORK2, tavg_O2_ZMIN,bid,k)

         call accumulate_tavg_field(WORK3, tavg_O2_ZMIN_DEPTH,bid,k)
      endif
   endif

   call accumulate_tavg_field(O2_PRODUCTION, tavg_O2_PRODUCTION,bid,k)

   call accumulate_tavg_field(O2_CONSUMPTION, tavg_O2_CONSUMPTION,bid,k)

   if (accumulate_tavg_now(tavg_AOU)) then
      WORK1 = O2SAT(TEMP, SALT, &
                    LAND_MASK(:,:,bid) .and. (k .le. KMT(:,:,bid)))
      WORK1 = WORK1 - p5*(TRACER_MODULE_OLD(:,:,k,o2_ind) + &
                          TRACER_MODULE_CUR(:,:,k,o2_ind))
      call accumulate_tavg_field(WORK1, tavg_AOU,bid,k)
   endif

   call accumulate_tavg_field(PAR_avg, tavg_PAR_avg,bid,k)

   call accumulate_tavg_field(graze_sp, tavg_graze_sp,bid,k)

   call accumulate_tavg_field(graze_diat, tavg_graze_diat,bid,k)

   call accumulate_tavg_field(graze_diaz, tavg_graze_diaz,bid,k)

   if (accumulate_tavg_now(tavg_graze_TOT)) then
      WORK1 = graze_sp + graze_diat + graze_diaz
      call accumulate_tavg_field(WORK1, tavg_graze_TOT,bid,k)
   endif

   call accumulate_tavg_field(sp_loss, tavg_sp_loss,bid,k)

   call accumulate_tavg_field(diat_loss, tavg_diat_loss,bid,k)

   call accumulate_tavg_field(diaz_loss, tavg_diaz_loss,bid,k)

   call accumulate_tavg_field(zoo_loss, tavg_zoo_loss,bid,k)

   call accumulate_tavg_field(sp_agg, tavg_sp_agg,bid,k)

   call accumulate_tavg_field(diat_agg, tavg_diat_agg,bid,k)

   call accumulate_tavg_field(photoC_sp, tavg_photoC_sp,bid,k)

   call accumulate_tavg_field(photoC_diat, tavg_photoC_diat,bid,k)

   call accumulate_tavg_field(photoC_diaz, tavg_photoC_diaz,bid,k)

   if (accumulate_tavg_now(tavg_photoC_TOT)) then
      WORK1 = photoC_sp + photoC_diat + photoC_diaz
      call accumulate_tavg_field(WORK1, tavg_photoC_TOT,bid,k)
   endif

   if (accumulate_tavg_now(tavg_photoC_sp_zint)) then
      if (partial_bottom_cells) then
         WORK1 = DZT(:,:,k,bid) * photoC_sp
      else
         WORK1 = dz(k) * photoC_sp
      endif
      call accumulate_tavg_field(WORK1, tavg_photoC_sp_zint,bid,k)
   endif

   if (accumulate_tavg_now(tavg_photoC_diat_zint)) then
      if (partial_bottom_cells) then
         WORK1 = DZT(:,:,k,bid) * photoC_diat
      else
         WORK1 = dz(k) * photoC_diat
      endif
      call accumulate_tavg_field(WORK1, tavg_photoC_diat_zint,bid,k)
   endif

   if (accumulate_tavg_now(tavg_photoC_diaz_zint)) then
      if (partial_bottom_cells) then
         WORK1 = DZT(:,:,k,bid) * photoC_diaz
      else
         WORK1 = dz(k) * photoC_diaz
      endif
      call accumulate_tavg_field(WORK1, tavg_photoC_diaz_zint,bid,k)
   endif

   if (accumulate_tavg_now(tavg_photoC_TOT_zint)) then
      if (partial_bottom_cells) then
         WORK1 = DZT(:,:,k,bid) * (photoC_sp + photoC_diat + photoC_diaz)
      else
         WORK1 = dz(k) * (photoC_sp + photoC_diat + photoC_diaz)
      endif
      call accumulate_tavg_field(WORK1, tavg_photoC_TOT_zint,bid,k)
   endif

   if (accumulate_tavg_now(tavg_photoC_NO3_sp) .or. &
       accumulate_tavg_now(tavg_photoC_NO3_sp_zint)) then
      where (VNtot_sp > c0)
         WORK1 = (VNO3_sp / VNtot_sp) * photoC_sp
      elsewhere
         WORK1 = c0
      end where

      call accumulate_tavg_field(WORK1, tavg_photoC_NO3_sp,bid,k)

      if (accumulate_tavg_now(tavg_photoC_NO3_sp_zint)) then
         if (partial_bottom_cells) then
            WORK1 = DZT(:,:,k,bid) * WORK1
         else
            WORK1 = dz(k) * WORK1
         endif
         call accumulate_tavg_field(WORK1, tavg_photoC_NO3_sp_zint,bid,k)
      endif
   endif

   if (accumulate_tavg_now(tavg_photoC_NO3_diat) .or. &
       accumulate_tavg_now(tavg_photoC_NO3_diat_zint)) then
      where (VNtot_diat > c0)
         WORK1 = (VNO3_diat / VNtot_diat) * photoC_diat
      elsewhere
         WORK1 = c0
      end where

      call accumulate_tavg_field(WORK1, tavg_photoC_NO3_diat,bid,k)

      if (accumulate_tavg_now(tavg_photoC_NO3_diat_zint)) then
         if (partial_bottom_cells) then
            WORK1 = DZT(:,:,k,bid) * WORK1
         else
            WORK1 = dz(k) * WORK1
         endif
         call accumulate_tavg_field(WORK1, tavg_photoC_NO3_diat_zint,bid,k)
      endif
   endif

   if (accumulate_tavg_now(tavg_photoC_NO3_diaz) .or. &
       accumulate_tavg_now(tavg_photoC_NO3_diaz_zint)) then
      where (VNtot_diaz > c0)
         WORK1 = (VNO3_diaz / VNtot_diaz) * photoC_diaz
      elsewhere
         WORK1 = c0
      end where

      call accumulate_tavg_field(WORK1, tavg_photoC_NO3_diaz,bid,k)

      if (accumulate_tavg_now(tavg_photoC_NO3_diaz_zint)) then
         if (partial_bottom_cells) then
            WORK1 = DZT(:,:,k,bid) * WORK1
         else
            WORK1 = dz(k) * WORK1
         endif
         call accumulate_tavg_field(WORK1, tavg_photoC_NO3_diaz_zint,bid,k)
      endif
   endif

   if (accumulate_tavg_now(tavg_photoC_NO3_TOT) .or. &
       accumulate_tavg_now(tavg_photoC_NO3_TOT_zint)) then
      where (VNtot_sp > c0)
         WORK1 = WORK1 + (VNO3_sp / VNtot_sp) * photoC_sp
      elsewhere
         WORK1 = c0
      end where
      where (VNtot_diat > c0)
         WORK1 = WORK1 + (VNO3_diat / VNtot_diat) * photoC_diat
      end where
      where (VNtot_diaz > c0)
         WORK1 = WORK1 + (VNO3_diaz / VNtot_diaz) * photoC_diaz
      end where

      call accumulate_tavg_field(WORK1, tavg_photoC_NO3_TOT,bid,k)

      if (accumulate_tavg_now(tavg_photoC_NO3_TOT_zint)) then
         if (partial_bottom_cells) then
            WORK1 = DZT(:,:,k,bid) * WORK1
         else
            WORK1 = dz(k) * WORK1
         endif
         call accumulate_tavg_field(WORK1, tavg_photoC_NO3_TOT_zint,bid,k)
      endif
   endif

   call accumulate_tavg_field(DOC_prod, tavg_DOC_prod,bid,k)

   call accumulate_tavg_field(DOC_remin, tavg_DOC_remin,bid,k)

   call accumulate_tavg_field(DON_prod, tavg_DON_prod,bid,k)

   call accumulate_tavg_field(DON_remin, tavg_DON_remin,bid,k)

   call accumulate_tavg_field(DOP_prod, tavg_DOP_prod,bid,k)

   call accumulate_tavg_field(DOP_remin, tavg_DOP_remin,bid,k)

   call accumulate_tavg_field(DOFe_prod, tavg_DOFe_prod,bid,k)

   call accumulate_tavg_field(DOFe_remin, tavg_DOFe_remin,bid,k)

   call accumulate_tavg_field(Fe_scavenge, tavg_Fe_scavenge,bid,k)

   call accumulate_tavg_field(Fe_scavenge_rate, tavg_Fe_scavenge_rate,bid,k)

   if (accumulate_tavg_now(tavg_CO3) .or. &
       accumulate_tavg_now(tavg_HCO3) .or. &
       accumulate_tavg_now(tavg_H2CO3) .or. &
       accumulate_tavg_now(tavg_pH_3D) .or. &
       accumulate_tavg_now(tavg_zsatcalc) .or. &
       accumulate_tavg_now(tavg_zsatarag)) then

      where (PH_PREV_3D(:,:,k,bid) /= c0)
         WORK1 = PH_PREV_3D(:,:,k,bid) - del_ph
         WORK2 = PH_PREV_3D(:,:,k,bid) + del_ph
      elsewhere
         WORK1 = phlo_3d_init
         WORK2 = phhi_3d_init
      end where

      call timer_start(ecosys_comp_CO3terms_timer, block_id=bid)
      call comp_CO3terms(bid, k, LAND_MASK(:,:,bid) .and. k <= KMT(:,:,bid), &
                         TEMP, SALT, DIC_loc, ALK_loc, PO4_loc, SiO3_loc, &
                         WORK1, WORK2, WORK3, H2CO3, HCO3, CO3)
      call timer_stop(ecosys_comp_CO3terms_timer, block_id=bid)
      call accumulate_tavg_field(CO3, tavg_CO3,bid,k)
      call accumulate_tavg_field(HCO3, tavg_HCO3,bid,k)
      call accumulate_tavg_field(H2CO3, tavg_H2CO3,bid,k)
      call accumulate_tavg_field(WORK3, tavg_pH_3D,bid,k)

      PH_PREV_3D(:,:,k,bid) = WORK3
   endif

   if (accumulate_tavg_now(tavg_co3_sat_calc) .or. &
       accumulate_tavg_now(tavg_zsatcalc) .or. &
       accumulate_tavg_now(tavg_co3_sat_arag) .or. &
       accumulate_tavg_now(tavg_zsatarag)) then
      call comp_co3_sat_vals(k, LAND_MASK(:,:,bid) .and. k <= KMT(:,:,bid), &
                             TEMP, SALT, WORK1, WORK2)
      call accumulate_tavg_field(WORK1, tavg_co3_sat_calc,bid,k)
      if (accumulate_tavg_now(tavg_zsatcalc)) then
         if (k == 1) then
            ! set to -1, i.e. depth not found yet,
            ! if mask == .true. and surface supersaturated to -1
            ZSATCALC(:,:,bid) = merge(-c1, c0, LAND_MASK(:,:,bid) .and. CO3 > WORK1)
         else
            where (ZSATCALC(:,:,bid) == -c1 .and. CO3 <= WORK1)
               ZSATCALC(:,:,bid) = zt(k-1) + (zt(k)-zt(k-1)) * &
                  CO3_CALC_ANOM_km1(:,:,bid) / (CO3_CALC_ANOM_km1(:,:,bid) - (CO3 - WORK1))
            endwhere
            where (ZSATCALC(:,:,bid) == -c1 .and. KMT(:,:,bid) == k)
               ZSATCALC(:,:,bid) = zw(k)
            endwhere
         endif
         CO3_CALC_ANOM_km1(:,:,bid) = CO3 - WORK1
         if (k == km) then
            call accumulate_tavg_field(ZSATCALC(:,:,bid), tavg_zsatcalc,bid,k)
         endif
      endif
      call accumulate_tavg_field(WORK2, tavg_co3_sat_arag,bid,k)
      if (accumulate_tavg_now(tavg_zsatarag)) then
         if (k == 1) then
            ! set to -1, i.e. depth not found yet,
            ! if mask == .true. and surface supersaturated to -1
            ZSATARAG(:,:,bid) = merge(-c1, c0, LAND_MASK(:,:,bid) .and. CO3 > WORK2)
         else
            where (ZSATARAG(:,:,bid) == -c1 .and. CO3 <= WORK2)
               ZSATARAG(:,:,bid) = zt(k-1) + (zt(k)-zt(k-1)) * &
                  CO3_ARAG_ANOM_km1(:,:,bid) / (CO3_ARAG_ANOM_km1(:,:,bid) - (CO3 - WORK2))
            endwhere
            where (ZSATARAG(:,:,bid) == -c1 .and. KMT(:,:,bid) == k)
               ZSATARAG(:,:,bid) = zw(k)
            endwhere
         endif
         CO3_ARAG_ANOM_km1(:,:,bid) = CO3 - WORK2
         if (k == km) then
            call accumulate_tavg_field(ZSATARAG(:,:,bid), tavg_zsatarag,bid,k)
         endif
      endif
   endif

   call timer_stop(ecosys_interior_timer, block_id=bid)

!-----------------------------------------------------------------------
!EOC

 end subroutine ecosys_set_interior

!***********************************************************************
!BOP
! !IROUTINE: init_particulate_terms
! !INTERFACE:

 subroutine init_particulate_terms(POC, P_CaCO3, P_SiO2, dust, P_iron, &
       QA_dust_def, NET_DUST_IN, this_block)

! !DESCRIPTION:
!  Set incoming fluxes (put into outgoing flux for first level usage).
!  Set dissolution length, production fraction and mass terms.
!
!  The first 6 arguments are intent(inout) in
!  order to preserve contents on other blocks.

! !INPUT/OUTPUT PARAMETERS:

   type(sinking_particle), intent(inout) :: &
      POC,          & ! base units = nmol C
      P_CaCO3,      & ! base units = nmol CaCO3
      P_SiO2,       & ! base units = nmol SiO2
      dust,         & ! base units = g
      P_iron          ! base units = nmol Fe

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic), intent(inout) :: &
      QA_dust_def     ! incoming deficit in the QA(dust) POC flux

! !INPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic), intent(in) :: &
      NET_DUST_IN     ! dust flux

   type (block), intent(in) :: &
      this_block      ! block info for the current block

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      bid                 ! local_block id

   character (char_len) :: &
      tracer_data_label   ! label for what is being updated

   character (char_len), dimension(:), allocatable :: &
      tracer_data_names   ! short names for input data fields

   integer (int_kind), dimension(:), allocatable :: &
      tracer_bndy_loc,  & ! location and field type for ghost
      tracer_bndy_type    !    cell updates

!-----------------------------------------------------------------------
!  parameters, from Armstrong et al. 2000
!
!  July 2002, length scale for excess POC and bSI modified by temperature
!  Value given here is at Tref of 30 deg. C, JKM
!-----------------------------------------------------------------------

    POC%diss      = 21000.0_r8  ! diss. length (cm), modified by TEMP
    POC%gamma     = c0          ! not used
    POC%mass      = 12.01_r8    ! molecular weight of POC
    POC%rho       = c0          ! not used

    P_CaCO3%diss  = 80000.0_r8  ! diss. length (cm)
    P_CaCO3%gamma = 0.55_r8     ! prod frac -> hard subclass
    P_CaCO3%mass  = 100.09_r8   ! molecular weight of CaCO
    P_CaCO3%rho   = 0.05_r8 * P_CaCO3%mass / POC%mass
                                       ! QA mass ratio for CaCO3
                                       ! This ratio is used in ecos_set_interior

    P_SiO2%diss   = 21000.0_r8  ! diss. length (cm), modified by TEMP
    P_SiO2%gamma  = 0.25_r8     ! prod frac -> hard subclass
    P_SiO2%mass   = 60.08_r8    ! molecular weight of SiO2
    P_SiO2%rho    = 0.05_r8 * P_SiO2%mass / POC%mass
                                       ! QA mass ratio for SiO2

    dust%diss     = 60000.0_r8  ! diss. length (cm)
    dust%gamma    = 0.97_r8     ! prod frac -> hard subclass
    dust%mass     = 1.0e9_r8    ! base units are already grams
    dust%rho      = 0.05_r8 * dust%mass / POC%mass
                                       ! QA mass ratio for dust

    P_iron%diss   = 60000.0_r8  ! diss. length (cm) - not used
    P_iron%gamma  = c0          ! prod frac -> hard subclass - not used
    P_iron%mass   = c0          ! not used
    P_iron%rho    = c0          ! not used

!-----------------------------------------------------------------------
!  Set incoming fluxes
!-----------------------------------------------------------------------

    bid = this_block%local_id

    P_CaCO3%sflux_out(:,:,bid) = c0
    P_CaCO3%hflux_out(:,:,bid) = c0

    P_SiO2%sflux_out(:,:,bid) = c0
    P_SiO2%hflux_out(:,:,bid) = c0

    if (dust_flux%has_data) then
       dust%sflux_out(:,:,bid) = (c1 - dust%gamma) * NET_DUST_IN(:,:,bid)
       dust%hflux_out(:,:,bid) = dust%gamma * NET_DUST_IN(:,:,bid)
    else
       dust%sflux_out(:,:,bid) = c0
       dust%hflux_out(:,:,bid) = c0
    endif

    P_iron%sflux_out(:,:,bid) = c0
    P_iron%hflux_out(:,:,bid) = c0

!-----------------------------------------------------------------------
!  Hard POC is QA flux and soft POC is excess POC.
!-----------------------------------------------------------------------

    POC%sflux_out(:,:,bid) = c0
    POC%hflux_out(:,:,bid) = c0

!-----------------------------------------------------------------------
!  Compute initial QA(dust) POC flux deficit.
!-----------------------------------------------------------------------

    QA_dust_def(:,:,bid) = dust%rho * &
       (dust%sflux_out(:,:,bid) + dust%hflux_out(:,:,bid))

!-----------------------------------------------------------------------
!EOC

 end subroutine init_particulate_terms

!***********************************************************************
!BOP
! !IROUTINE: compute_particulate_terms
! !INTERFACE:

 subroutine compute_particulate_terms(k, POC, P_CaCO3, P_SiO2, dust, P_iron, &
       QA_dust_def, TEMP, O2_loc, this_block)

! !DESCRIPTION:
!  Compute outgoing fluxes and remineralization terms. Assumes that
!  production terms have been set. Incoming fluxes are assumed to be the
!  outgoing fluxes from the previous level.
!
!  It is assumed that there is no production of dust.
!
!  Instantaneous remineralization in the bottom cell is implemented by
!  setting the outgoing flux to zero.
!
!  For POC, the hard subclass is the POC flux qualitatively associated
!  with the ballast flux. The soft subclass is the excess POC flux.
!
!  Remineralization for the non-iron particulate pools is computing
!  by first computing the outgoing flux and then computing the
!  remineralization from conservation, i.e.
!     flux_in - flux_out + prod * dz - remin * dz == 0.
!
!  For iron, remineralization is first computed from POC remineralization
!  and then flux_out is computed from conservation. If the resulting
!  flux_out is negative or should be zero because of the sea floor, the
!  remineralization is adjusted.
!  Note: all the sinking iron is in the P_iron%sflux pool, hflux Fe not
!        explicitly tracked, it is assumed that total iron remin is
!        proportional to total POC remin.
!
!  Based upon Armstrong et al. 2000
!
!  July 2002, added temperature effect on remin length scale of
!  excess POC (all soft POM& Iron) and on SiO2.
!  new variable passed into ballast, Tfunc, main Temperature function
!  computed in ecosystem routine.  scaling factor for dissolution
!  of excess POC, Fe, and Bsi now varies with location (f(temperature)).
!
!  Added diffusive iron flux from sediments at depths < 1100m,
!  based on Johnson et al., 1999, value of 5 umolFe/m2/day,
!      this value too high, using 2 umolFe/m2/day here
!
!  Allow hard fraction of ballast to remin with long length scale 40,000m
!     thus ~ 10% of hard ballast remins over 4000m water column.
!
!  Sinking dust flux is decreased by assumed instant solubility/dissolution
!     at ocean surface from the parm_Fe_bioavail.
!
!  Modified to allow different Q10 factors for soft POM and bSI remin,
!  water TEMP is now passed in instead of Tfunc (1/2005, JKM)

! !USES:

#ifdef CCSMCOUPLED
   use shr_sys_mod, only: shr_sys_abort
#endif

! !INPUT PARAMETERS:

   integer (int_kind), intent(in) :: k ! vertical model level

   real (r8), dimension(nx_block,ny_block), intent(in) :: &
      TEMP,         & ! temperature for scaling functions
      O2_loc          ! dissolved oxygen used to change/increase POC%diss

   type (block), intent(in) :: &
      this_block      ! block info for the current block

! !INPUT/OUTPUT PARAMETERS:

   type(sinking_particle), intent(inout) :: &
      POC,          & ! base units = nmol C
      P_CaCO3,      & ! base units = nmol CaCO3
      P_SiO2,       & ! base units = nmol SiO2
      dust,         & ! base units = g
      P_iron          ! base units = nmol Fe

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic), intent(inout) :: &
      QA_dust_def     ! incoming deficit in the QA(dust) POC flux

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   real (r8) :: poc_diss ! diss. length used (cm)

   character(*), parameter :: &
      subname = 'ecosys_mod:compute_particulate_terms'

   real (r8), dimension(nx_block,ny_block) :: &
      WORK,               & ! temporary for summed quantities to be averaged
      TfuncP,             & ! temperature scaling from soft POM remin
      TfuncS,             & ! temperature scaling from soft POM remin
      DECAY_CaCO3,        & ! scaling factor for dissolution of CaCO3
      DECAY_dust,         & ! scaling factor for dissolution of dust
      DECAY_Hard,         & ! scaling factor for dissolution of Hard Ballast
      DECAY_HardDust        ! scaling factor for dissolution of Hard dust

   real (r8) :: &
      decay_POC_E,        & ! scaling factor for dissolution of excess POC
      decay_SiO2,         & ! scaling factor for dissolution of SiO2
      POC_PROD_avail,     & ! POC production available for excess POC flux
      new_QA_dust_def,    & ! outgoing deficit in the QA(dust) POC flux
      dz_loc, dzr_loc       ! dz, dzr at a particular i,j location

   integer (int_kind) :: &
      i, j,               & ! loop indices
      bid                   ! local_block id

   logical (log_kind) :: &
      poc_error             ! POC error flag

!-----------------------------------------------------------------------
!  incoming fluxes are outgoing fluxes from previous level
!-----------------------------------------------------------------------

   bid = this_block%local_id

   P_CaCO3%sflux_in(:,:,bid) = P_CaCO3%sflux_out(:,:,bid)
   P_CaCO3%hflux_in(:,:,bid) = P_CaCO3%hflux_out(:,:,bid)

   P_SiO2%sflux_in(:,:,bid) = P_SiO2%sflux_out(:,:,bid)
   P_SiO2%hflux_in(:,:,bid) = P_SiO2%hflux_out(:,:,bid)

   dust%sflux_in(:,:,bid) = dust%sflux_out(:,:,bid)
   dust%hflux_in(:,:,bid) = dust%hflux_out(:,:,bid)

   POC%sflux_in(:,:,bid) = POC%sflux_out(:,:,bid)
   POC%hflux_in(:,:,bid) = POC%hflux_out(:,:,bid)

   P_iron%sflux_in(:,:,bid) = P_iron%sflux_out(:,:,bid)
   P_iron%hflux_in(:,:,bid) = P_iron%hflux_out(:,:,bid)

!-----------------------------------------------------------------------
!  compute decay factors
!-----------------------------------------------------------------------

   if (partial_bottom_cells) then
      DECAY_CaCO3    = exp(-DZT(:,:,k,bid) / P_CaCO3%diss)
      DECAY_dust     = exp(-DZT(:,:,k,bid) / dust%diss)
      DECAY_Hard     = exp(-DZT(:,:,k,bid) / 4.0e6_r8)
      DECAY_HardDust = exp(-DZT(:,:,k,bid) / 1.2e7_r8)
   else
      DECAY_CaCO3    = exp(-dz(k) / P_CaCO3%diss)
      DECAY_dust     = exp(-dz(k) / dust%diss)
      DECAY_Hard     = exp(-dz(k) / 4.0e6_r8)
      DECAY_HardDust = exp(-dz(k) / 1.2e7_r8)
   endif

!----------------------------------------------------------------------
!   Tref = 30.0 reference temperature (deg. C)
!
!   Using q10 formulation with Q10 value of 1.1 soft POM (TfuncP) and
!       a Q10 value of 2.5 soft bSi (TfuncS)
!-----------------------------------------------------------------------

!
!  NOTE: Turning off temperature effect on POM lengthscale, three instances
!    of TfuncP below have been removed, see comment lines.
!-------------------------------------------------------------------------
!  TfuncP = 1.1_r8**(((TEMP + T0_Kelvin) - (Tref + T0_Kelvin)) / c10)

   TfuncS = 2.5_r8**(((TEMP + T0_Kelvin) - (Tref + T0_Kelvin)) / c10)

   poc_error = .false.

   do j = 1,ny_block
      do i = 1,nx_block

         if (LAND_MASK(i,j,bid) .and. k <= KMT(i,j,bid)) then

            if (partial_bottom_cells) then
               dz_loc = DZT(i,j,k,bid)
            else
               dz_loc = dz(k)
            endif
            dzr_loc = c1 / dz_loc

!-----------------------------------------------------------------------
!  increase POC diss where O2 concentrations are low
!-----------------------------------------------------------------------

!            if (O2_loc(i,j) > c0) then
!               poc_diss = min((POC%diss * 2.0_r8), (POC%diss / &
!                  (O2_loc(i,j) / (O2_loc(i,j) + 5.0_r8))))
!            else
!               poc_diss = POC%diss * 2.0_r8
!            endif

             poc_diss = POC%diss

!-----------------------------------------------------------------------
!  decay_POC_E and decay_SiO2 set locally, modified by Tfunc
!-----------------------------------------------------------------------

!            decay_POC_E = exp(-dz_loc / (poc_diss / TfuncP(i,j)))

            decay_POC_E = exp(-dz_loc / poc_diss)
            decay_SiO2  = exp(-dz_loc / (P_SiO2%diss / TfuncS(i,j)))

!-----------------------------------------------------------------------
!  Set outgoing fluxes for non-iron pools.
!  The outoing fluxes for ballast materials are from the
!  solution of the coresponding continuous ODE across the model
!  level. The ODE has a constant source term and linear decay.
!  It is assumed that there is no sub-surface dust production.
!-----------------------------------------------------------------------

            P_CaCO3%sflux_out(i,j,bid) = P_CaCO3%sflux_in(i,j,bid) * DECAY_CaCO3(i,j) + &
               P_CaCO3%prod(i,j,bid) * ((c1 - P_CaCO3%gamma) * (c1 - DECAY_CaCO3(i,j)) &
                  * P_CaCO3%diss)

            P_CaCO3%hflux_out(i,j,bid) = P_CaCO3%hflux_in(i,j,bid) * DECAY_Hard(i,j) + &
               P_CaCO3%prod(i,j,bid) * (P_CaCO3%gamma * dz_loc)

            P_SiO2%sflux_out(i,j,bid) = P_SiO2%sflux_in(i,j,bid) * decay_SiO2 + &
               P_SiO2%prod(i,j,bid) * ((c1 - P_SiO2%gamma) * (c1 - decay_SiO2) &
                  * (P_SiO2%diss / TfuncS(i,j)))

            P_SiO2%hflux_out(i,j,bid) = P_SiO2%hflux_in(i,j,bid) * DECAY_Hard(i,j) + &
               P_SiO2%prod(i,j,bid) * (P_SiO2%gamma * dz_loc)

            dust%sflux_out(i,j,bid) = dust%sflux_in(i,j,bid) * DECAY_dust(i,j)

            dust%hflux_out(i,j,bid) = dust%hflux_in(i,j,bid) * DECAY_HardDust(i,j)

!-----------------------------------------------------------------------
!  Compute how much POC_PROD is available for deficit reduction
!  and excess POC flux after subtracting off fraction of non-dust
!  ballast production from net POC_PROD.
!-----------------------------------------------------------------------

            POC_PROD_avail = POC%prod(i,j,bid) - &
               P_CaCO3%rho * P_CaCO3%prod(i,j,bid) - &
               P_SiO2%rho * P_SiO2%prod(i,j,bid)

!-----------------------------------------------------------------------
!  Check for POC production bounds violations
!-----------------------------------------------------------------------

            if (POC_PROD_avail < c0) then
               poc_error = .true.
            endif

!-----------------------------------------------------------------------
!  Compute 1st approximation to new QA_dust_def, the QA_dust
!  deficit leaving the cell. Ignore POC_PROD_avail at this stage.
!-----------------------------------------------------------------------

            if (QA_dust_def(i,j,bid) > 0) then
               new_QA_dust_def = QA_dust_def(i,j,bid) * &
                  (dust%sflux_out(i,j,bid) + dust%hflux_out(i,j,bid)) / &
                  (dust%sflux_in(i,j,bid) + dust%hflux_in(i,j,bid))
            else
               new_QA_dust_def = c0
            endif

!-----------------------------------------------------------------------
!  Use POC_PROD_avail to reduce new_QA_dust_def.
!-----------------------------------------------------------------------

            if (new_QA_dust_def > c0) then
               new_QA_dust_def = new_QA_dust_def - POC_PROD_avail * dz_loc
               if (new_QA_dust_def < c0) then
                  POC_PROD_avail = -new_QA_dust_def * dzr_loc
                  new_QA_dust_def = c0
               else
                  POC_PROD_avail = c0
               endif
            endif

            QA_dust_def(i,j,bid) = new_QA_dust_def

!-----------------------------------------------------------------------
!  Compute outgoing POC fluxes. QA POC flux is computing using
!  ballast fluxes and new_QA_dust_def. If no QA POC flux came in
!  and no production occured, then no QA POC flux goes out. This
!  shortcut is present to avoid roundoff cancellation errors from
!  the dust%rho * dust_flux_out - QA_dust_def computation.
!  Any POC_PROD_avail still remaining goes into excess POC flux.
!-----------------------------------------------------------------------

            if (POC%hflux_in(i,j,bid) == c0 .and. POC%prod(i,j,bid) == c0) then
               POC%hflux_out(i,j,bid) = c0
            else
               POC%hflux_out(i,j,bid) = P_CaCO3%rho * &
                  (P_CaCO3%sflux_out(i,j,bid) + P_CaCO3%hflux_out(i,j,bid)) + &
                  P_SiO2%rho * &
                  (P_SiO2%sflux_out(i,j,bid) + P_SiO2%hflux_out(i,j,bid)) + &
                  dust%rho * &
                  (dust%sflux_out(i,j,bid) + dust%hflux_out(i,j,bid)) - &
                  new_QA_dust_def
               POC%hflux_out(i,j,bid) = max(POC%hflux_out(i,j,bid), c0)
            endif

            POC%sflux_out(i,j,bid) = POC%sflux_in(i,j,bid) * decay_POC_E + &
               POC_PROD_avail *((c1 - decay_POC_E) * &
               poc_diss)

!               (poc_diss / TfuncP(i,j)))

!-----------------------------------------------------------------------
!  Compute remineralization terms. It is assumed that there is no
!  sub-surface dust production.
!-----------------------------------------------------------------------

            P_CaCO3%remin(i,j,bid) = P_CaCO3%prod(i,j,bid) + &
               ((P_CaCO3%sflux_in(i,j,bid) - P_CaCO3%sflux_out(i,j,bid)) + &
               (P_CaCO3%hflux_in(i,j,bid) - P_CaCO3%hflux_out(i,j,bid))) * dzr_loc

            P_SiO2%remin(i,j,bid) = P_SiO2%prod(i,j,bid) + &
               ((P_SiO2%sflux_in(i,j,bid) - P_SiO2%sflux_out(i,j,bid)) + &
               (P_SiO2%hflux_in(i,j,bid) - P_SiO2%hflux_out(i,j,bid))) * dzr_loc

            POC%remin(i,j,bid) = POC%prod(i,j,bid) + &
               ((POC%sflux_in(i,j,bid) - POC%sflux_out(i,j,bid)) + &
               (POC%hflux_in(i,j,bid) - POC%hflux_out(i,j,bid))) * dzr_loc

            dust%remin(i,j,bid) = &
               ((dust%sflux_in(i,j,bid) - dust%sflux_out(i,j,bid)) + &
               (dust%hflux_in(i,j,bid) - dust%hflux_out(i,j,bid))) * dzr_loc

!-----------------------------------------------------------------------
!  Compute iron remineralization and flux out.
!-----------------------------------------------------------------------

            if (POC%sflux_in(i,j,bid) + POC%hflux_in(i,j,bid) == c0) then
               P_iron%remin(i,j,bid) = (POC%remin(i,j,bid) * parm_Red_Fe_C)
            else
               P_iron%remin(i,j,bid) = (POC%remin(i,j,bid) * &
                  (P_iron%sflux_in(i,j,bid) + P_iron%hflux_in(i,j,bid)) / &
                  (POC%sflux_in(i,j,bid) + POC%hflux_in(i,j,bid))) + &
                  (P_iron%sflux_in(i,j,bid) * 3.0e-6_r8)
            endif

            P_iron%sflux_out(i,j,bid) = P_iron%sflux_in(i,j,bid) + dz_loc * &
               ((c1 - P_iron%gamma) * P_iron%prod(i,j,bid) - P_iron%remin(i,j,bid))

            if (P_iron%sflux_out(i,j,bid) < c0) then
               P_iron%sflux_out(i,j,bid) = c0
               P_iron%remin(i,j,bid) = P_iron%sflux_in(i,j,bid) * dzr_loc + &
                  (c1 - P_iron%gamma) * P_iron%prod(i,j,bid)
            endif

!-----------------------------------------------------------------------
!  Compute iron release from dust remin/dissolution
!
!  dust remin gDust = 0.035 / 55.847 * 1.0e9 = 626712.0 nmolFe
!                      gFe     molFe     nmolFe
!  Also add in Fe source from sediments if applicable to this cell.
!-----------------------------------------------------------------------


            P_iron%remin(i,j,bid) = P_iron%remin(i,j,bid) &
               + dust%remin(i,j,bid) * dust_to_Fe &
               + (FESEDFLUX(i,j,k,bid) * dzr_loc)

!maltrud what is up here--jkm does not have this
            P_iron%hflux_out(i,j,bid) = P_iron%hflux_in(i,j,bid)

         else
            P_CaCO3%sflux_out(i,j,bid) = c0
            P_CaCO3%hflux_out(i,j,bid) = c0
            P_CaCO3%remin(i,j,bid) = c0

            P_SiO2%sflux_out(i,j,bid) = c0
            P_SiO2%hflux_out(i,j,bid) = c0
            P_SiO2%remin(i,j,bid) = c0

            dust%sflux_out(i,j,bid) = c0
            dust%hflux_out(i,j,bid) = c0
            dust%remin(i,j,bid) = c0

            POC%sflux_out(i,j,bid) = c0
            POC%hflux_out(i,j,bid) = c0
            POC%remin(i,j,bid) = c0

            P_iron%sflux_out(i,j,bid) = c0
            P_iron%hflux_out(i,j,bid) = c0
            P_iron%remin(i,j,bid) = c0
         endif

!-----------------------------------------------------------------------
!  Remineralize everything in bottom cell.
!-----------------------------------------------------------------------

         if (LAND_MASK(i,j,bid) .and. k == KMT(i,j,bid)) then
            P_CaCO3%remin(i,j,bid) = P_CaCO3%remin(i,j,bid) + &
               (P_CaCO3%sflux_out(i,j,bid) + P_CaCO3%hflux_out(i,j,bid)) * dzr_loc
            P_CaCO3%sflux_out(i,j,bid) = c0
            P_CaCO3%hflux_out(i,j,bid) = c0

            P_SiO2%remin(i,j,bid) = P_SiO2%remin(i,j,bid) + &
               (P_SiO2%sflux_out(i,j,bid) + P_SiO2%hflux_out(i,j,bid)) * dzr_loc
            P_SiO2%sflux_out(i,j,bid) = c0
            P_SiO2%hflux_out(i,j,bid) = c0

!           dust%remin(i,j,bid) = dust%remin(i,j,bid) + &
!              (dust%sflux_out(i,j,bid) + dust%hflux_out(i,j,bid)) * dzr_loc
!           dust%sflux_out(i,j,bid) = c0
!           dust%hflux_out(i,j,bid) = c0

            POC%remin(i,j,bid) = POC%remin(i,j,bid) + &
               (POC%sflux_out(i,j,bid) + POC%hflux_out(i,j,bid)) * dzr_loc
            POC%sflux_out(i,j,bid) = c0
            POC%hflux_out(i,j,bid) = c0

            P_iron%remin(i,j,bid) = P_iron%remin(i,j,bid) + &
               (P_iron%sflux_out(i,j,bid) + P_iron%hflux_out(i,j,bid)) * dzr_loc
            P_iron%sflux_out(i,j,bid) = c0
            P_iron%hflux_out(i,j,bid) = c0

         endif

      end do
   end do

#ifdef CCSMCOUPLED
    if (poc_error) then
       call shr_sys_abort(subname /&
                       &/ ': mass ratio of ballast ' /&
                       &/ 'production exceeds POC production')
    endif
#endif

!-----------------------------------------------------------------------
!  Set tavg variables.
!-----------------------------------------------------------------------

   if (accumulate_tavg_now(tavg_POC_FLUX_IN)) then
      WORK = POC%sflux_in(:,:,bid) + POC%hflux_in(:,:,bid)
      call accumulate_tavg_field(WORK, tavg_POC_FLUX_IN,bid,k)
   endif

   call accumulate_tavg_field(POC%prod(:,:,bid), tavg_POC_PROD,bid,k)

   call accumulate_tavg_field(POC%remin(:,:,bid), tavg_POC_REMIN,bid,k)

   if (accumulate_tavg_now(tavg_CaCO3_FLUX_IN)) then
      WORK = P_CaCO3%sflux_in(:,:,bid) + P_CaCO3%hflux_in(:,:,bid)
      call accumulate_tavg_field(WORK, tavg_CaCO3_FLUX_IN,bid,k)
   endif

   call accumulate_tavg_field(P_CaCO3%prod(:,:,bid), tavg_CaCO3_PROD,bid,k)

   call accumulate_tavg_field(P_CaCO3%remin(:,:,bid), tavg_CaCO3_REMIN,bid,k)

   if (accumulate_tavg_now(tavg_SiO2_FLUX_IN)) then
      WORK = P_SiO2%sflux_in(:,:,bid) + P_SiO2%hflux_in(:,:,bid)
      call accumulate_tavg_field(WORK, tavg_SiO2_FLUX_IN,bid,k)
   endif

   call accumulate_tavg_field(P_SiO2%prod(:,:,bid), tavg_SiO2_PROD,bid,k)

   call accumulate_tavg_field(P_SiO2%remin(:,:,bid), tavg_SiO2_REMIN,bid,k)

   if (accumulate_tavg_now(tavg_dust_FLUX_IN)) then
      WORK = dust%sflux_in(:,:,bid) + dust%hflux_in(:,:,bid)
      call accumulate_tavg_field(WORK, tavg_dust_FLUX_IN,bid,k)
   endif

   call accumulate_tavg_field(dust%remin(:,:,bid), tavg_dust_REMIN,bid,k)

   if (accumulate_tavg_now(tavg_P_iron_FLUX_IN)) then
      WORK = P_iron%sflux_in(:,:,bid) + P_iron%hflux_in(:,:,bid)
      call accumulate_tavg_field(WORK, tavg_P_iron_FLUX_IN,bid,k)
   endif

   call accumulate_tavg_field(P_iron%prod(:,:,bid), tavg_P_iron_PROD,bid,k)

   call accumulate_tavg_field(P_iron%remin(:,:,bid), tavg_P_iron_REMIN,bid,k)

!-----------------------------------------------------------------------
!EOC

 end subroutine compute_particulate_terms

!***********************************************************************
!BOP
! !IROUTINE: ecosys_init_sflux
! !INTERFACE:

 subroutine ecosys_init_sflux

! !DESCRIPTION:
!  Initialize surface flux computations for ecosys tracer module.
!
! !REVISION HISTORY:
!  same as module

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   character(*), parameter :: subname = 'ecosys_mod:ecosys_init_sflux'

   logical (log_kind) :: &
      luse_INTERP_WORK     ! does INTERP_WORK need to be allocated

   integer (int_kind) :: &
      n,                 & ! index for looping over tracers
      iblock               ! index for looping over blocks

   real (r8), dimension (nx_block,ny_block) :: WORK

   real (r8), dimension (nx_block,ny_block,12,max_blocks_clinic), target :: &
      WORK_READ            ! temporary space to read in fields

!-----------------------------------------------------------------------

   luse_INTERP_WORK = .false.

!-----------------------------------------------------------------------
!  read gas flux forcing (if required)
!-----------------------------------------------------------------------

   if ((lflux_gas_o2 .or. lflux_gas_co2) .and. &
       gas_flux_forcing_iopt == gas_flux_forcing_iopt_file) then

      luse_INTERP_WORK = .true.

!-----------------------------------------------------------------------
!  first, read ice file
!-----------------------------------------------------------------------

      allocate(fice_file%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
      if (trim(fice_file%input%filename) == 'unknown') &
         fice_file%input%filename = gas_flux_forcing_file

      call read_field(fice_file%input%file_fmt, &
                      fice_file%input%filename, &
                      fice_file%input%file_varname, &
                      WORK_READ)

      !$OMP PARALLEL DO PRIVATE(iblock, n)
      do iblock=1,nblocks_clinic
      do n=1,12
         fice_file%DATA(:,:,iblock,1,n) = WORK_READ(:,:,n,iblock)
         where (.not. LAND_MASK(:,:,iblock)) &
            fice_file%DATA(:,:,iblock,1,n) = c0
         fice_file%DATA(:,:,iblock,1,n) = &
            fice_file%DATA(:,:,iblock,1,n) * fice_file%input%scale_factor
      end do
      end do
      !$OMP END PARALLEL DO

      call find_forcing_times(fice_file%data_time, &
                              fice_file%data_inc, fice_file%interp_type, &
                              fice_file%data_next, fice_file%data_time_min_loc, &
                              fice_file%data_update, fice_file%data_type)

!-----------------------------------------------------------------------
!  next, read piston velocity file
!-----------------------------------------------------------------------

      allocate(xkw_file%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
      if (trim(xkw_file%input%filename) == 'unknown') &
         xkw_file%input%filename = gas_flux_forcing_file

      call read_field(xkw_file%input%file_fmt, &
                      xkw_file%input%filename, &
                      xkw_file%input%file_varname, &
                      WORK_READ)

      !$OMP PARALLEL DO PRIVATE(iblock, n)
      do iblock=1,nblocks_clinic
      do n=1,12
         xkw_file%DATA(:,:,iblock,1,n) = WORK_READ(:,:,n,iblock)
         where (.not. LAND_MASK(:,:,iblock)) &
            xkw_file%DATA(:,:,iblock,1,n) = c0
         xkw_file%DATA(:,:,iblock,1,n) = &
            xkw_file%DATA(:,:,iblock,1,n) * xkw_file%input%scale_factor
      end do
      end do
      !$OMP END PARALLEL DO

      call find_forcing_times(xkw_file%data_time, &
                              xkw_file%data_inc, xkw_file%interp_type, &
                              xkw_file%data_next, xkw_file%data_time_min_loc, &
                              xkw_file%data_update, xkw_file%data_type)

!-----------------------------------------------------------------------
!  last, read atmospheric pressure file
!-----------------------------------------------------------------------

      allocate(ap_file%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
      if (trim(ap_file%input%filename) == 'unknown') &
         ap_file%input%filename = gas_flux_forcing_file

      call read_field(ap_file%input%file_fmt, &
                      ap_file%input%filename, &
                      ap_file%input%file_varname, &
                      WORK_READ)

      !$OMP PARALLEL DO PRIVATE(iblock, n)
      do iblock=1,nblocks_clinic
      do n=1,12
         ap_file%DATA(:,:,iblock,1,n) = WORK_READ(:,:,n,iblock)
         where (.not. LAND_MASK(:,:,iblock)) &
            ap_file%DATA(:,:,iblock,1,n) = c0
         ap_file%DATA(:,:,iblock,1,n) = &
            ap_file%DATA(:,:,iblock,1,n) * ap_file%input%scale_factor
      end do
      end do
      !$OMP END PARALLEL DO

      call find_forcing_times(ap_file%data_time, &
                              ap_file%data_inc, ap_file%interp_type, &
                              ap_file%data_next, ap_file%data_time_min_loc, &
                              ap_file%data_update, ap_file%data_type)

    endif

!-----------------------------------------------------------------------
!  load dust flux fields (if required)
!-----------------------------------------------------------------------

   if (trim(dust_flux%input%filename) /= 'none' .and. &
       trim(dust_flux%input%filename) /= 'unknown') then

      luse_INTERP_WORK = .true.

      allocate(dust_flux%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
      if (trim(dust_flux%input%filename) == 'unknown') &
         dust_flux%input%filename = gas_flux_forcing_file

      call read_field(dust_flux%input%file_fmt, &
                      dust_flux%input%filename, &
                      dust_flux%input%file_varname, &
                      WORK_READ)

      !$OMP PARALLEL DO PRIVATE(iblock, n)
      do iblock=1,nblocks_clinic
      do n=1,12
         dust_flux%DATA(:,:,iblock,1,n) = WORK_READ(:,:,n,iblock)
         where (.not. LAND_MASK(:,:,iblock)) &
            dust_flux%DATA(:,:,iblock,1,n) = c0
         dust_flux%DATA(:,:,iblock,1,n) = &
            dust_flux%DATA(:,:,iblock,1,n) * dust_flux%input%scale_factor
      end do
      end do
      !$OMP END PARALLEL DO

      call find_forcing_times(dust_flux%data_time, &
                              dust_flux%data_inc, dust_flux%interp_type, &
                              dust_flux%data_next, dust_flux%data_time_min_loc, &
                              dust_flux%data_update, dust_flux%data_type)

      dust_flux%has_data = .true.
   else
      dust_flux%has_data = .false.
   endif

!-----------------------------------------------------------------------
!  load iron flux fields (if required)
!-----------------------------------------------------------------------

   if (trim(iron_flux%input%filename) /= 'none' .and. &
       trim(iron_flux%input%filename) /= 'unknown') then

      luse_INTERP_WORK = .true.

      allocate(iron_flux%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
      if (trim(iron_flux%input%filename) == 'unknown') &
         iron_flux%input%filename = gas_flux_forcing_file

      call read_field(iron_flux%input%file_fmt, &
                      iron_flux%input%filename, &
                      iron_flux%input%file_varname, &
                      WORK_READ)

      !$OMP PARALLEL DO PRIVATE(iblock, n)
      do iblock=1,nblocks_clinic
      do n=1,12
         iron_flux%DATA(:,:,iblock,1,n) = WORK_READ(:,:,n,iblock)
         where (.not. LAND_MASK(:,:,iblock)) &
            iron_flux%DATA(:,:,iblock,1,n) = c0
         iron_flux%DATA(:,:,iblock,1,n) = &
            iron_flux%DATA(:,:,iblock,1,n) * iron_flux%input%scale_factor
      end do
      end do
      !$OMP END PARALLEL DO

      call find_forcing_times(iron_flux%data_time, &
                              iron_flux%data_inc, iron_flux%interp_type, &
                              iron_flux%data_next, iron_flux%data_time_min_loc, &
                              iron_flux%data_update, iron_flux%data_type)

      iron_flux%has_data = .true.
   else
      iron_flux%has_data = .false.
   endif

!-----------------------------------------------------------------------
!  load nox & noy flux fields (if required)
!-----------------------------------------------------------------------

   if (trim(ndep_data_type) /= 'none' .and. &
       trim(ndep_data_type) /= 'monthly-calendar' .and. &
       trim(ndep_data_type) /= 'shr_stream') then
      call document(subname, 'ndep_data_type', ndep_data_type)
      call exit_POP(sigAbort, 'unknown ndep_data_type')
   endif

   if (trim(ndep_data_type) == 'monthly-calendar' .and. &
       trim(nox_flux_monthly%input%filename) /= 'none' .and. &
       trim(nox_flux_monthly%input%filename) /= 'unknown') then

      luse_INTERP_WORK = .true.

      allocate(nox_flux_monthly%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
      if (trim(nox_flux_monthly%input%filename) == 'unknown') &
         nox_flux_monthly%input%filename = gas_flux_forcing_file

      call read_field(nox_flux_monthly%input%file_fmt, &
                      nox_flux_monthly%input%filename, &
                      nox_flux_monthly%input%file_varname, &
                      WORK_READ)

      !$OMP PARALLEL DO PRIVATE(iblock, n)
      do iblock=1,nblocks_clinic
      do n=1,12
         nox_flux_monthly%DATA(:,:,iblock,1,n) = WORK_READ(:,:,n,iblock)
         where (.not. LAND_MASK(:,:,iblock)) &
            nox_flux_monthly%DATA(:,:,iblock,1,n) = c0
         nox_flux_monthly%DATA(:,:,iblock,1,n) = &
            nox_flux_monthly%DATA(:,:,iblock,1,n) * nox_flux_monthly%input%scale_factor
      end do
      end do
      !$OMP END PARALLEL DO

      call find_forcing_times(nox_flux_monthly%data_time, &
                              nox_flux_monthly%data_inc, nox_flux_monthly%interp_type, &
                              nox_flux_monthly%data_next, nox_flux_monthly%data_time_min_loc, &
                              nox_flux_monthly%data_update, nox_flux_monthly%data_type)

      nox_flux_monthly%has_data = .true.
   else
      nox_flux_monthly%has_data = .false.
   endif

   if (trim(ndep_data_type) == 'monthly-calendar' .and. &
       trim(nhy_flux_monthly%input%filename) /= 'none' .and. &
       trim(nhy_flux_monthly%input%filename) /= 'unknown') then

      luse_INTERP_WORK = .true.

      allocate(nhy_flux_monthly%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
      if (trim(nhy_flux_monthly%input%filename) == 'unknown') &
         nhy_flux_monthly%input%filename = gas_flux_forcing_file

      call read_field(nhy_flux_monthly%input%file_fmt, &
                      nhy_flux_monthly%input%filename, &
                      nhy_flux_monthly%input%file_varname, &
                      WORK_READ)

      !$OMP PARALLEL DO PRIVATE(iblock, n)
      do iblock=1,nblocks_clinic
      do n=1,12
         nhy_flux_monthly%DATA(:,:,iblock,1,n) = WORK_READ(:,:,n,iblock)
         where (.not. LAND_MASK(:,:,iblock)) &
            nhy_flux_monthly%DATA(:,:,iblock,1,n) = c0
         nhy_flux_monthly%DATA(:,:,iblock,1,n) = &
            nhy_flux_monthly%DATA(:,:,iblock,1,n) * nhy_flux_monthly%input%scale_factor
      end do
      end do
      !$OMP END PARALLEL DO

      call find_forcing_times(nhy_flux_monthly%data_time, &
                              nhy_flux_monthly%data_inc, nhy_flux_monthly%interp_type, &
                              nhy_flux_monthly%data_next, nhy_flux_monthly%data_time_min_loc, &
                              nhy_flux_monthly%data_update, nhy_flux_monthly%data_type)

      nhy_flux_monthly%has_data = .true.
   else
      nhy_flux_monthly%has_data = .false.
   endif

!-----------------------------------------------------------------------

#ifndef CCSMCOUPLED
   if (trim(ndep_data_type) == 'shr_stream') then
      call document(subname, 'ndep_data_type', ndep_data_type)
      call exit_POP(sigAbort, &
         'shr_stream option only supported when CCSMCOUPLED is defined')
   endif
#endif

!-----------------------------------------------------------------------
!  allocate space for interpolate_forcing
!-----------------------------------------------------------------------

   if (luse_INTERP_WORK) &
      allocate(INTERP_WORK(nx_block,ny_block,max_blocks_clinic,1))

!-----------------------------------------------------------------------
!  load iron PATCH flux fields (if required)
!-----------------------------------------------------------------------

   if (liron_patch) then

!maltrud iron patch
!  assume patch file has same normalization and format as deposition file

      allocate(IRON_PATCH_FLUX(nx_block,ny_block,max_blocks_clinic))

      if (trim(iron_flux%input%filename) == 'unknown') &
         iron_flux%input%filename = gas_flux_forcing_file

      call read_field(iron_flux%input%file_fmt, &
                      iron_flux%input%filename, &
                      iron_patch_flux_filename, &
                      IRON_PATCH_FLUX)

      !$OMP PARALLEL DO PRIVATE(iblock, n)
      do iblock=1,nblocks_clinic
      do n=1,12
         where (.not. LAND_MASK(:,:,iblock)) &
            IRON_PATCH_FLUX(:,:,iblock) = c0
         iron_flux%DATA(:,:,iblock,1,n) = &
            IRON_PATCH_FLUX(:,:,iblock) * iron_flux%input%scale_factor
      end do
      end do
      !$OMP END PARALLEL DO

   endif

!-----------------------------------------------------------------------
!  register and set
!     fco2, the air-sea co2 gas flux
!-----------------------------------------------------------------------

   call named_field_register('SFLUX_CO2', sflux_co2_nf_ind)
   !$OMP PARALLEL DO PRIVATE(iblock,WORK)
   do iblock=1,nblocks_clinic
      WORK = c0
      call named_field_set(sflux_co2_nf_ind, iblock, WORK)
   end do
   !$OMP END PARALLEL DO

!-----------------------------------------------------------------------
!  verify running coupled if gas fluxes use coupler forcing
!-----------------------------------------------------------------------

   if ((lflux_gas_o2 .or. lflux_gas_co2) .and. &
       (gas_flux_forcing_iopt == gas_flux_forcing_iopt_drv .or. &
        atm_co2_iopt == atm_co2_iopt_drv_prog .or. &
        atm_co2_iopt == atm_co2_iopt_drv_diag) .and. &
       .not. registry_match('lcoupled')) then
      call exit_POP(sigAbort, 'ecosys_init: ecosys module requires the ' /&
                           &/ 'flux coupler when gas_flux_forcing_opt=drv')
   endif

!-----------------------------------------------------------------------
!  get named field index for atmospheric CO2, if required
!-----------------------------------------------------------------------

   if (lflux_gas_co2 .and. atm_co2_iopt == atm_co2_iopt_drv_prog) then
      call named_field_get_index('ATM_CO2_PROG', atm_co2_nf_ind, &
                                 exit_on_err=.false.)
      if (atm_co2_nf_ind == 0) then
         call exit_POP(sigAbort, 'ecosys_init: ecosys module requires ' /&
                              &/ 'atmopsheric CO2 from the flux coupler ' /&
                              &/ 'and it is not present')
      endif
   endif

   if (lflux_gas_co2 .and. atm_co2_iopt == atm_co2_iopt_drv_diag) then
      call named_field_get_index('ATM_CO2_DIAG', atm_co2_nf_ind, &
                                 exit_on_err=.false.)
      if (atm_co2_nf_ind == 0) then
         call exit_POP(sigAbort, 'ecosys_init: ecosys module requires ' /&
                              &/ 'atmopsheric CO2 from the flux coupler ' /&
                              &/ 'and it is not present')
      endif
   endif

!-----------------------------------------------------------------------
!EOC

 end subroutine ecosys_init_sflux

!***********************************************************************
!BOP
! !IROUTINE: ecosys_init_interior_restore
! !INTERFACE:

 subroutine ecosys_init_interior_restore

! !DESCRIPTION:
!  Initialize interior restoring computations for ecosys tracer module.
!
! !REVISION HISTORY:
!  same as module

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      n,                   & ! index for looping over tracers
      k,                   & ! index for looping over levels
      i,j,                 & ! index for looping over horiz. dims.
      iblock                 ! index for looping over blocks

   real (r8) :: &
      subsurf_fesed          ! sum of subsurface fesed values

!-----------------------------------------------------------------------
!  initialize restoring timescale (if required)
!-----------------------------------------------------------------------

   if (lrest_po4 .or. lrest_no3 .or. lrest_sio3) then

      if (lnutr_variable_restore) then

         allocate(NUTR_RESTORE_RTAU(nx_block,ny_block,max_blocks_clinic))
         allocate(NUTR_RESTORE_MAX_LEVEL(nx_block,ny_block,max_blocks_clinic))

         call read_field(nutr_variable_rest_file_fmt, &
                         nutr_variable_rest_file,     &
                         'NUTR_RESTORE_MAX_LEVEL',    &
                         NUTR_RESTORE_RTAU)

         NUTR_RESTORE_MAX_LEVEL = nint(NUTR_RESTORE_RTAU)

         call read_field(nutr_variable_rest_file_fmt, &
                         nutr_variable_rest_file,     &
                         'NUTR_RESTORE_RTAU',         &
                         NUTR_RESTORE_RTAU)

         NUTR_RESTORE_RTAU = NUTR_RESTORE_RTAU/seconds_in_day ! convert days to secs

      else

         do k = 1,km
            if (zt(k) < rest_z0) then
               nutr_rest_time_inv(k) = rest_time_inv_surf
            else if (zt(k) > rest_z1) then
               nutr_rest_time_inv(k) = rest_time_inv_deep
            else if (rest_z1 == rest_z0) then
               nutr_rest_time_inv(k) = rest_time_inv_surf + p5 * &
                    (rest_time_inv_deep - rest_time_inv_surf)
            else
               nutr_rest_time_inv(k) = rest_time_inv_surf + &
                    (zt(k) - rest_z0) / (rest_z1 - rest_z0) * &
                    (rest_time_inv_deep - rest_time_inv_surf)
            endif
         end do

      endif  !  variable restoring

   endif  ! restoring

!-----------------------------------------------------------------------
!  load restoring fields (if required)
!-----------------------------------------------------------------------

   if (lrest_po4) then

      allocate(PO4_CLIM(nx_block,ny_block,km,max_blocks_clinic))

      call read_field(po4_rest%file_fmt, &
                      po4_rest%filename, &
                      po4_rest%file_varname, &
                      PO4_CLIM)

      do iblock=1,nblocks_clinic
         do k = 1, km
            where (.not. LAND_MASK(:,:,iblock) .or. k > KMT(:,:,iblock)) &
               PO4_CLIM(:,:,k,iblock) = c0
            PO4_CLIM(:,:,k,iblock) = PO4_CLIM(:,:,k,iblock) * po4_rest%scale_factor
         enddo
      end do

   endif

   if (lrest_no3) then

      allocate(NO3_CLIM(nx_block,ny_block,km,max_blocks_clinic))

      call read_field(no3_rest%file_fmt, &
                      no3_rest%filename, &
                      no3_rest%file_varname, &
                      NO3_CLIM)

      do iblock=1,nblocks_clinic
         do k = 1, km
            where (.not. LAND_MASK(:,:,iblock) .or. k > KMT(:,:,iblock)) &
               NO3_CLIM(:,:,k,iblock) = c0
            NO3_CLIM(:,:,k,iblock) = NO3_CLIM(:,:,k,iblock) * no3_rest%scale_factor
         enddo
      end do

   endif

   if (lrest_sio3) then

      allocate(SiO3_CLIM(nx_block,ny_block,km,max_blocks_clinic))

      call read_field(sio3_rest%file_fmt, &
                      sio3_rest%filename, &
                      sio3_rest%file_varname, &
                      SiO3_CLIM)

      do iblock=1,nblocks_clinic
         do k = 1, km
            where (.not. LAND_MASK(:,:,iblock) .or. k > KMT(:,:,iblock)) &
               SiO3_CLIM(:,:,k,iblock) = c0
            SiO3_CLIM(:,:,k,iblock) = SiO3_CLIM(:,:,k,iblock) * sio3_rest%scale_factor
         enddo
      end do

   endif

!-----------------------------------------------------------------------
!  load fesedflux
!  add subsurface positives to 1 level shallower, to accomodate overflow pop-ups
!-----------------------------------------------------------------------

   allocate(FESEDFLUX(nx_block,ny_block,km,max_blocks_clinic))

   call read_field(fesedflux_input%file_fmt, &
                   fesedflux_input%filename, &
                   fesedflux_input%file_varname, &
                   FESEDFLUX)

   do iblock=1,nblocks_clinic
      do j=1,ny_block
      do i=1,nx_block
         if (KMT(i,j,iblock) > 0 .and. KMT(i,j,iblock) < km) then
            subsurf_fesed = c0
            do k=KMT(i,j,iblock)+1,km
               subsurf_fesed = subsurf_fesed + FESEDFLUX(i,j,k,iblock)
            enddo
            FESEDFLUX(i,j,KMT(i,j,iblock),iblock) = FESEDFLUX(i,j,KMT(i,j,iblock),iblock) + subsurf_fesed
         endif
      enddo
      enddo

      do k = 1, km
         where (.not. LAND_MASK(:,:,iblock) .or. k > KMT(:,:,iblock)) &
            FESEDFLUX(:,:,k,iblock) = c0
         FESEDFLUX(:,:,k,iblock) = FESEDFLUX(:,:,k,iblock) * fesedflux_input%scale_factor
      enddo
   end do

!-----------------------------------------------------------------------
!EOC

 end subroutine ecosys_init_interior_restore

!***********************************************************************
!BOP
! !IROUTINE: ecosys_set_sflux
! !INTERFACE:

 subroutine ecosys_set_sflux(SHF_QSW_RAW, SHF_QSW, &
                             U10_SQR,IFRAC,PRESS,SST,SSS, &
                             SURF_VALS_OLD,SURF_VALS_CUR,STF_MODULE)
   use seq_io_mod, only : seq_io_getiotype, seq_io_getiosys

! !DESCRIPTION:
!  Compute surface fluxes for ecosys tracer module.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic), intent(in) :: &
      SHF_QSW_RAW,  &! penetrative solar heat flux, from coupler (degC*cm/s)
      SHF_QSW,      &! SHF_QSW used by physics, may have diurnal cylce imposed (degC*cm/s)
      U10_SQR,      &! 10m wind speed squared (cm/s)**2
      IFRAC,        &! sea ice fraction (non-dimensional)
      PRESS,        &! sea level atmospheric pressure (dyne/cm**2)
      SST,          &! sea surface temperature (C)
      SSS            ! sea surface salinity (psu)

   real (r8), dimension(nx_block,ny_block,ecosys_tracer_cnt,max_blocks_clinic), &
      intent(in) :: SURF_VALS_OLD, SURF_VALS_CUR ! module tracers

! !INPUT/OUTPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block,ecosys_tracer_cnt,max_blocks_clinic), &
      intent(inout) :: STF_MODULE

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   character(*), parameter :: subname = 'ecosys_mod:ecosys_set_sflux'

   logical (log_kind) :: first_call = .true.

   type (block) :: &
      this_block      ! block info for the current block

   integer (int_kind) :: &
      i,j,iblock,n, & ! loop indices
      mcdate,sec,   & ! date vals for shr_strdata_advance
      errorCode       ! errorCode from HaloUpdate call

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
      IFRAC_USED,   & ! used ice fraction (non-dimensional)
      XKW_USED,     & ! portion of piston velocity (cm/s)
      AP_USED,      & ! used atm pressure (converted from dyne/cm**2 to atm)
      IRON_FLUX_IN    ! iron flux

   real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
      SHR_STREAM_WORK

   real (r8), dimension(nx_block,ny_block) :: &
      XKW_ICE,      & ! common portion of piston vel., (1-fice)*xkw (cm/s)
      SCHMIDT_USED, & ! used Schmidt number
      PV,           & ! piston velocity (cm/s)
      O2SAT_1atm,   & ! O2 saturation @ 1 atm (mmol/m^3)
      O2SAT_USED,   & ! used O2 saturation (mmol/m^3)
      XCO2,         & ! atmospheric co2 conc. (dry-air, 1 atm)
      FLUX            ! tracer flux (nmol/cm^2/s)

   real (r8), dimension(nx_block) :: &
      PHLO,         & ! lower bound for ph in solver
      PHHI,         & ! upper bound for ph in solver
      PH_NEW,       & ! computed PH from solver
      DIC_ROW,      & ! row of DIC values for solver
      ALK_ROW,      & ! row of ALK values for solver
      PO4_ROW,      & ! row of PO4 values for solver
      SiO3_ROW,     & ! row of SiO3 values for solver
      CO2STAR_ROW,  & ! CO2STAR from solver
      DCO2STAR_ROW, & ! DCO2STAR from solver
      pCO2SURF_ROW, & ! pCO2SURF from solver
      DpCO2_ROW       ! DpCO2 from solver

   character (char_len) :: &
      tracer_data_label,       & ! label for what is being updated
      ndep_shr_stream_fldList

   character (char_len), dimension(1) :: &
      tracer_data_names          ! short names for input data fields

   integer (int_kind), dimension(1) :: &
      tracer_bndy_loc,         & ! location and field type for ghost
      tracer_bndy_type           !    cell updates

   real (r8), dimension(nx_block,ny_block) :: &
      WORK1, WORK2 ! temporaries for averages

   real (r8) :: scalar_temp

!-----------------------------------------------------------------------
!  local parameters
!-----------------------------------------------------------------------

   real (r8), parameter :: &
      xkw_coeff = 8.6e-9_r8     ! a = 0.31 cm/hr s^2/m^2 in (s/cm)

!-----------------------------------------------------------------------

   call timer_start(ecosys_sflux_timer)

!-----------------------------------------------------------------------

   if (check_time_flag(comp_surf_avg_flag))  &
      call comp_surf_avg(SURF_VALS_OLD,SURF_VALS_CUR)

!-----------------------------------------------------------------------
!  fluxes initially set to 0
!  set Chl field for short-wave absorption
!  store incoming shortwave in PAR_out field, converting to W/m^2
!-----------------------------------------------------------------------

   scalar_temp = f_qsw_par / hflux_factor

   !$OMP PARALLEL DO PRIVATE(iblock,WORK1)
   do iblock = 1, nblocks_clinic
      STF_MODULE(:,:,:,iblock) = c0

      WORK1 = max(c0,p5*(SURF_VALS_OLD(:,:,spChl_ind,iblock) + &
                         SURF_VALS_CUR(:,:,spChl_ind,iblock))) + &
              max(c0,p5*(SURF_VALS_OLD(:,:,diatChl_ind,iblock) + &
                         SURF_VALS_CUR(:,:,diatChl_ind,iblock))) + &
              max(c0,p5*(SURF_VALS_OLD(:,:,diazChl_ind,iblock) + &
                         SURF_VALS_CUR(:,:,diazChl_ind,iblock)))
      call named_field_set(totChl_surf_nf_ind, iblock, WORK1)

      if (ecosys_qsw_distrb_const) then
         PAR_out(:,:,iblock) = SHF_QSW_RAW(:,:,iblock)
      else
         PAR_out(:,:,iblock) = SHF_QSW(:,:,iblock)
      endif

      where (LAND_MASK(:,:,iblock))
         PAR_out(:,:,iblock) = max(c0, scalar_temp * PAR_out(:,:,iblock))
      elsewhere
         PAR_out(:,:,iblock) = c0
      end where
   enddo
   !$OMP END PARALLEL DO

!-----------------------------------------------------------------------
!  Interpolate gas flux forcing data if necessary
!-----------------------------------------------------------------------

   if ((lflux_gas_o2 .or. lflux_gas_co2) .and. &
        gas_flux_forcing_iopt == gas_flux_forcing_iopt_file) then
      if (thour00 >= fice_file%data_update) then
         tracer_data_names = fice_file%input%file_varname
         tracer_bndy_loc   = field_loc_center
         tracer_bndy_type  = field_type_scalar
         tracer_data_label = 'Ice Fraction'
         call update_forcing_data(fice_file%data_time,      &
            fice_file%data_time_min_loc,  fice_file%interp_type,    &
            fice_file%data_next,          fice_file%data_update,    &
            fice_file%data_type,          fice_file%data_inc,       &
            fice_file%DATA(:,:,:,:,1:12), fice_file%data_renorm,    &
            tracer_data_label,            tracer_data_names,        &
            tracer_bndy_loc,              tracer_bndy_type,         &
            fice_file%filename,           fice_file%input%file_fmt)
      endif
      call interpolate_forcing(INTERP_WORK, &
         fice_file%DATA(:,:,:,:,1:12), &
         fice_file%data_time,         fice_file%interp_type, &
         fice_file%data_time_min_loc, fice_file%interp_freq, &
         fice_file%interp_inc,        fice_file%interp_next, &
         fice_file%interp_last,       0)
      IFRAC_USED = INTERP_WORK(:,:,:,1)

      if (thour00 >= xkw_file%data_update) then
         tracer_data_names = xkw_file%input%file_varname
         tracer_bndy_loc   = field_loc_center
         tracer_bndy_type  = field_type_scalar
         tracer_data_label = 'Piston Velocity'
         call update_forcing_data(xkw_file%data_time,      &
            xkw_file%data_time_min_loc,  xkw_file%interp_type,    &
            xkw_file%data_next,          xkw_file%data_update,    &
            xkw_file%data_type,          xkw_file%data_inc,       &
            xkw_file%DATA(:,:,:,:,1:12), xkw_file%data_renorm,    &
            tracer_data_label,           tracer_data_names,       &
            tracer_bndy_loc,             tracer_bndy_type,        &
            xkw_file%filename,           xkw_file%input%file_fmt)
      endif
      call interpolate_forcing(INTERP_WORK,     &
         xkw_file%DATA(:,:,:,:,1:12), &
         xkw_file%data_time,         xkw_file%interp_type, &
         xkw_file%data_time_min_loc, xkw_file%interp_freq, &
         xkw_file%interp_inc,        xkw_file%interp_next, &
         xkw_file%interp_last,       0)
      XKW_USED = INTERP_WORK(:,:,:,1)

      if (thour00 >= ap_file%data_update) then
         tracer_data_names = ap_file%input%file_varname
         tracer_bndy_loc   = field_loc_center
         tracer_bndy_type  = field_type_scalar
         tracer_data_label = 'Atmospheric Pressure'
         call update_forcing_data(ap_file%data_time,    &
            ap_file%data_time_min_loc,  ap_file%interp_type,    &
            ap_file%data_next,          ap_file%data_update,    &
            ap_file%data_type,          ap_file%data_inc,       &
            ap_file%DATA(:,:,:,:,1:12), ap_file%data_renorm,    &
            tracer_data_label,          tracer_data_names,      &
            tracer_bndy_loc,            tracer_bndy_type,       &
            ap_file%filename,           ap_file%input%file_fmt)
      endif
      call interpolate_forcing(INTERP_WORK, &
         ap_file%DATA(:,:,:,:,1:12), &
         ap_file%data_time,         ap_file%interp_type, &
         ap_file%data_time_min_loc, ap_file%interp_freq, &
         ap_file%interp_inc,        ap_file%interp_next, &
         ap_file%interp_last,       0)
      AP_USED = INTERP_WORK(:,:,:,1)

   endif

!-----------------------------------------------------------------------
!  calculate gas flux quantities if necessary
!-----------------------------------------------------------------------

   if (lflux_gas_o2 .or. lflux_gas_co2) then

      !$OMP PARALLEL DO PRIVATE(iblock,j,XKW_ICE,SCHMIDT_USED,PV,O2SAT_USED, &
      !$OMP                     O2SAT_1atm,FLUX,XCO2,PHLO,PHHI,DIC_ROW,ALK_ROW, &
      !$OMP                     PO4_ROW,SiO3_ROW,PH_NEW,CO2STAR_ROW, &
      !$OMP                     DCO2STAR_ROW,pCO2SURF_ROW,DpCO2_ROW)

      do iblock = 1, nblocks_clinic

!-----------------------------------------------------------------------
!  Apply OCMIP ice fraction mask when input is from a file.
!-----------------------------------------------------------------------

         if (gas_flux_forcing_iopt == gas_flux_forcing_iopt_file) then
            where (IFRAC_USED(:,:,iblock) < 0.2000_r8) &
               IFRAC_USED(:,:,iblock) = 0.2000_r8
            where (IFRAC_USED(:,:,iblock) > 0.9999_r8) &
               IFRAC_USED(:,:,iblock) = 0.9999_r8
         endif

         if (gas_flux_forcing_iopt == gas_flux_forcing_iopt_drv) then
            IFRAC_USED(:,:,iblock) = IFRAC(:,:,iblock)
            where (IFRAC_USED(:,:,iblock) < c0) IFRAC_USED(:,:,iblock) = c0
            where (IFRAC_USED(:,:,iblock) > c1) IFRAC_USED(:,:,iblock) = c1
            XKW_USED(:,:,iblock) = xkw_coeff * U10_SQR(:,:,iblock)
            AP_USED(:,:,iblock) = PRESS(:,:,iblock)
         endif

!-----------------------------------------------------------------------
!  assume PRESS is in cgs units (dyne/cm**2) since that is what is
!    required for pressure forcing in barotropic
!  want units to be atmospheres
!  convertion from dyne/cm**2 to Pascals is P(mks) = P(cgs)/10.
!  convertion from Pascals to atm is P(atm) = P(Pa)/101.325e+3_r8
!-----------------------------------------------------------------------

         AP_USED(:,:,iblock) = PRESS(:,:,iblock) / 101.325e+4_r8

         ECO_SFLUX_TAVG(:,:,1,iblock) = IFRAC_USED(:,:,iblock)
         ECO_SFLUX_TAVG(:,:,2,iblock) = XKW_USED(:,:,iblock)
         ECO_SFLUX_TAVG(:,:,3,iblock) = AP_USED(:,:,iblock)

!-----------------------------------------------------------------------
!  Compute XKW_ICE. XKW is zero over land, so XKW_ICE is too.
!-----------------------------------------------------------------------

         XKW_ICE = (c1 - IFRAC_USED(:,:,iblock)) * XKW_USED(:,:,iblock)

!-----------------------------------------------------------------------
!  compute O2 flux
!-----------------------------------------------------------------------

         if (lflux_gas_o2) then
            SCHMIDT_USED = SCHMIDT_O2(SST(:,:,iblock), LAND_MASK(:,:,iblock))

            O2SAT_1atm = O2SAT(SST(:,:,iblock),SSS(:,:,iblock),   &
                               LAND_MASK(:,:,iblock))

            where (LAND_MASK(:,:,iblock))
               PV = XKW_ICE * SQRT(660.0_r8 / SCHMIDT_USED)
               O2SAT_USED = AP_USED(:,:,iblock) * O2SAT_1atm
               FLUX = PV * (O2SAT_USED - p5*(SURF_VALS_OLD(:,:,o2_ind,iblock) + &
                                             SURF_VALS_CUR(:,:,o2_ind,iblock)))
               STF_MODULE(:,:,o2_ind,iblock) = FLUX
            elsewhere
               O2SAT_USED = c0
               STF_MODULE(:,:,o2_ind,iblock) = c0
            end where

            ECO_SFLUX_TAVG(:,:,4,iblock) = PV
            ECO_SFLUX_TAVG(:,:,5,iblock) = SCHMIDT_USED
            ECO_SFLUX_TAVG(:,:,6,iblock) = O2SAT_USED
            ECO_SFLUX_TAVG(:,:,16,iblock) = STF_MODULE(:,:,o2_ind,iblock)

         endif  ! lflux_gas_o2

!-----------------------------------------------------------------------
!  compute CO2 flux, computing disequilibrium one row at a time
!-----------------------------------------------------------------------

         if (lflux_gas_co2) then

            SCHMIDT_USED = SCHMIDT_CO2(SST(:,:,iblock), LAND_MASK(:,:,iblock))

            where (LAND_MASK(:,:,iblock))
               PV = XKW_ICE * SQRT(660.0_r8 / SCHMIDT_USED)
            elsewhere
               PV = c0
            end where

!-----------------------------------------------------------------------
!  Set XCO2
!-----------------------------------------------------------------------

            select case (atm_co2_iopt)
            case (atm_co2_iopt_const)
               XCO2 = atm_co2_const
            case (atm_co2_iopt_drv_prog, atm_co2_iopt_drv_diag)
               call named_field_get(atm_co2_nf_ind, iblock, XCO2)
            end select

            do j = 1,ny_block
               where (PH_PREV(:,j,iblock) /= c0)
                  PHLO = PH_PREV(:,j,iblock) - del_ph
                  PHHI = PH_PREV(:,j,iblock) + del_ph
               elsewhere
                  PHLO = phlo_surf_init
                  PHHI = phhi_surf_init
               end where

               DIC_ROW = p5*(SURF_VALS_OLD(:,j,dic_ind,iblock) + &
                             SURF_VALS_CUR(:,j,dic_ind,iblock))
               ALK_ROW = p5*(SURF_VALS_OLD(:,j,alk_ind,iblock) + &
                             SURF_VALS_CUR(:,j,alk_ind,iblock))
               PO4_ROW = p5*(SURF_VALS_OLD(:,j,po4_ind,iblock) + &
                             SURF_VALS_CUR(:,j,po4_ind,iblock))
               SiO3_ROW = p5*(SURF_VALS_OLD(:,j,sio3_ind,iblock) + &
                              SURF_VALS_CUR(:,j,sio3_ind,iblock))

               call co2calc_row(iblock, j, LAND_MASK(:,j,iblock), locmip_k1_k2_bug_fix, &
                                SST(:,j,iblock), SSS(:,j,iblock), &
                                DIC_ROW, ALK_ROW, PO4_ROW, SiO3_ROW, &
                                PHLO, PHHI, PH_NEW, XCO2(:,j), &
                                AP_USED(:,j,iblock), CO2STAR_ROW, &
                                DCO2STAR_ROW, pCO2SURF_ROW, DpCO2_ROW)

               PH_PREV(:,j,iblock) = PH_NEW

               FLUX(:,j) = PV(:,j) * DCO2STAR_ROW

               ECO_SFLUX_TAVG(:,j, 7,iblock) = CO2STAR_ROW
               ECO_SFLUX_TAVG(:,j, 8,iblock) = DCO2STAR_ROW
               ECO_SFLUX_TAVG(:,j, 9,iblock) = pCO2SURF_ROW
               ECO_SFLUX_TAVG(:,j,10,iblock) = DpCO2_ROW

            end do

!-----------------------------------------------------------------------
!  set air-sea co2 gas flux named field, converting units from
!  nmol/cm^2/s (positive down) to kg CO2/m^2/s (positive down)
!-----------------------------------------------------------------------

            call named_field_set(sflux_co2_nf_ind, iblock, 44.0e-8_r8 * FLUX)

            STF_MODULE(:,:,dic_ind,iblock) = STF_MODULE(:,:,dic_ind,iblock) + FLUX

            ECO_SFLUX_TAVG(:,:,11,iblock) = PV
            ECO_SFLUX_TAVG(:,:,12,iblock) = SCHMIDT_USED
            ECO_SFLUX_TAVG(:,:,13,iblock) = FLUX
            ECO_SFLUX_TAVG(:,:,14,iblock) = PH_PREV(:,:,iblock)
            ECO_SFLUX_TAVG(:,:,15,iblock) = XCO2

         endif  !  lflux_gas_co2

      enddo
      !$OMP END PARALLEL DO

   endif  ! lflux_gas_o2 .or. lflux_gas_co2

!-----------------------------------------------------------------------
!  calculate iron and dust fluxes if necessary
!-----------------------------------------------------------------------

   if (iron_flux%has_data) then
      if (thour00 >= iron_flux%data_update) then
         tracer_data_names = iron_flux%input%file_varname
         tracer_bndy_loc   = field_loc_center
         tracer_bndy_type  = field_type_scalar
         tracer_data_label = 'Iron Flux'
         call update_forcing_data(iron_flux%data_time,    &
            iron_flux%data_time_min_loc,  iron_flux%interp_type,    &
            iron_flux%data_next,          iron_flux%data_update,    &
            iron_flux%data_type,          iron_flux%data_inc,       &
            iron_flux%DATA(:,:,:,:,1:12), iron_flux%data_renorm,    &
            tracer_data_label,            tracer_data_names,        &
            tracer_bndy_loc,              tracer_bndy_type,         &
            iron_flux%filename,           iron_flux%input%file_fmt)
      endif
      call interpolate_forcing(INTERP_WORK,     &
         iron_flux%DATA(:,:,:,:,1:12), &
         iron_flux%data_time,         iron_flux%interp_type, &
         iron_flux%data_time_min_loc, iron_flux%interp_freq, &
         iron_flux%interp_inc,        iron_flux%interp_next, &
         iron_flux%interp_last,       0)
      if (liron_patch .and. imonth == iron_patch_month) then
         IRON_FLUX_IN = INTERP_WORK(:,:,:,1) + IRON_PATCH_FLUX
      else
         IRON_FLUX_IN = INTERP_WORK(:,:,:,1)
      endif
   else
      IRON_FLUX_IN = c0
   endif

   IRON_FLUX_IN = IRON_FLUX_IN * parm_Fe_bioavail

   STF_MODULE(:,:,fe_ind,:) = IRON_FLUX_IN

   if (dust_flux%has_data) then
      if (thour00 >= dust_flux%data_update) then
         tracer_data_names = dust_flux%input%file_varname
         tracer_bndy_loc   = field_loc_center
         tracer_bndy_type  = field_type_scalar
         tracer_data_label = 'Dust Flux'
         call update_forcing_data(dust_flux%data_time,    &
            dust_flux%data_time_min_loc,  dust_flux%interp_type,    &
            dust_flux%data_next,          dust_flux%data_update,    &
            dust_flux%data_type,          dust_flux%data_inc,       &
            dust_flux%DATA(:,:,:,:,1:12), dust_flux%data_renorm,    &
            tracer_data_label,            tracer_data_names,        &
            tracer_bndy_loc,              tracer_bndy_type,         &
            dust_flux%filename,           dust_flux%input%file_fmt)
      endif
      call interpolate_forcing(INTERP_WORK, &
         dust_flux%DATA(:,:,:,:,1:12),    &
         dust_flux%data_time,         dust_flux%interp_type, &
         dust_flux%data_time_min_loc, dust_flux%interp_freq, &
         dust_flux%interp_inc,        dust_flux%interp_next, &
         dust_flux%interp_last,       0)
      dust_FLUX_IN = INTERP_WORK(:,:,:,1)

!-----------------------------------------------------------------------
!  Reduce surface dust flux due to assumed instant surface dissolution
!-----------------------------------------------------------------------

      dust_FLUX_IN = dust_FLUX_IN * (c1 - parm_Fe_bioavail)
   else
      dust_FLUX_IN = c0
   endif

!-----------------------------------------------------------------------
!  calculate nox and nhy fluxes if necessary
!-----------------------------------------------------------------------

   if (nox_flux_monthly%has_data) then
      if (thour00 >= nox_flux_monthly%data_update) then
         tracer_data_names = nox_flux_monthly%input%file_varname
         tracer_bndy_loc   = field_loc_center
         tracer_bndy_type  = field_type_scalar
         tracer_data_label = 'NOx Flux'
         call update_forcing_data(nox_flux_monthly%data_time,    &
            nox_flux_monthly%data_time_min_loc,  nox_flux_monthly%interp_type, &
            nox_flux_monthly%data_next,          nox_flux_monthly%data_update, &
            nox_flux_monthly%data_type,          nox_flux_monthly%data_inc,    &
            nox_flux_monthly%DATA(:,:,:,:,1:12), nox_flux_monthly%data_renorm, &
            tracer_data_label,                   tracer_data_names,            &
            tracer_bndy_loc,                     tracer_bndy_type,             &
            nox_flux_monthly%filename,           nox_flux_monthly%input%file_fmt)
      endif
      call interpolate_forcing(INTERP_WORK,     &
         nox_flux_monthly%DATA(:,:,:,:,1:12), &
         nox_flux_monthly%data_time,         nox_flux_monthly%interp_type, &
         nox_flux_monthly%data_time_min_loc, nox_flux_monthly%interp_freq, &
         nox_flux_monthly%interp_inc,        nox_flux_monthly%interp_next, &
         nox_flux_monthly%interp_last,       0)
      STF_MODULE(:,:,no3_ind,:) = INTERP_WORK(:,:,:,1)
   endif

   if (nhy_flux_monthly%has_data) then
      if (thour00 >= nhy_flux_monthly%data_update) then
         tracer_data_names = nhy_flux_monthly%input%file_varname
         tracer_bndy_loc   = field_loc_center
         tracer_bndy_type  = field_type_scalar
         tracer_data_label = 'NHy Flux'
         call update_forcing_data(nhy_flux_monthly%data_time,    &
            nhy_flux_monthly%data_time_min_loc,  nhy_flux_monthly%interp_type, &
            nhy_flux_monthly%data_next,          nhy_flux_monthly%data_update, &
            nhy_flux_monthly%data_type,          nhy_flux_monthly%data_inc,    &
            nhy_flux_monthly%DATA(:,:,:,:,1:12), nhy_flux_monthly%data_renorm, &
            tracer_data_label,                   tracer_data_names,            &
            tracer_bndy_loc,                     tracer_bndy_type,             &
            nhy_flux_monthly%filename,           nhy_flux_monthly%input%file_fmt)
      endif
      call interpolate_forcing(INTERP_WORK,     &
         nhy_flux_monthly%DATA(:,:,:,:,1:12), &
         nhy_flux_monthly%data_time,         nhy_flux_monthly%interp_type, &
         nhy_flux_monthly%data_time_min_loc, nhy_flux_monthly%interp_freq, &
         nhy_flux_monthly%interp_inc,        nhy_flux_monthly%interp_next, &
         nhy_flux_monthly%interp_last,       0)
      STF_MODULE(:,:,nh4_ind,:) = INTERP_WORK(:,:,:,1)
   endif

#ifdef CCSMCOUPLED
   if (trim(ndep_data_type) == 'shr_stream') then
      if (first_call) then

         ndep_shr_stream_fldList = ' '
         do n = 1, ndep_shr_stream_var_cnt
            if (n == ndep_shr_stream_no_ind) &
               ndep_shr_stream_fldList = trim(ndep_shr_stream_fldList) /&
                  &/ 'NOy_deposition'
            if (n == ndep_shr_stream_nh_ind) &
               ndep_shr_stream_fldList = trim(ndep_shr_stream_fldList) /&
                  &/ 'NHx_deposition'
            if (n < ndep_shr_stream_var_cnt) &
               ndep_shr_stream_fldList = trim(ndep_shr_stream_fldList) /&
                  &/ ':'
         end do

         call shr_strdata_create(ndep_sdat,name='ndep data',                   &
                                 mpicom=POP_communicator,                      &
                                 compid=POP_MCT_OCNID,                         &
                                 gsmap=POP_MCT_gsMap_o, ggrid=POP_MCT_dom_o,   &
                                 nxg=nx_global, nyg=ny_global,                 &
                                 yearFirst=ndep_shr_stream_year_first,         &
                                 yearLast=ndep_shr_stream_year_last,           &
                                 yearAlign=ndep_shr_stream_year_align,         &
                                 offset=0,                                     &
                                 domFilePath='',                               &
                                 domFileName=ndep_shr_stream_file,             &
                                 domTvarName='time',                           &
                                 domXvarName='TLONG', domYvarName='TLAT',      &
                                 domAreaName='TAREA', domMaskName='KMT',       &
                                 FilePath='',                                  &
                                 FileName=(/trim(ndep_shr_stream_file)/),      &
                                 fldListFile=ndep_shr_stream_fldList,          &
                                 fldListModel=ndep_shr_stream_fldList,         &
                                 pio_subsystem=seq_io_getiosys('OCN'),         &
                                 pio_iotype=seq_io_getiotype('OCN'),           &
                                 fillalgo='none', mapalgo='none')
         if (my_task == master_task) then
            call shr_strdata_print(ndep_sdat)
         endif
         first_call = .false.
      endif

      mcdate = iyear*10000 + imonth*100 + iday
      sec = isecond + 60 * (iminute + 60 * ihour)

      call timer_start(ecosys_shr_strdata_advance_timer)
      call shr_strdata_advance(ndep_sdat, mcdate, sec, &
                               POP_communicator, 'ndep')
      call timer_stop(ecosys_shr_strdata_advance_timer)

      !
      ! process NO3 flux, store results in SHR_STREAM_WORK array
      ! instead of directly into STF_MODULE
      ! to avoid argument copies in HaloUpdate calls
      !
      n = 0
      do iblock = 1, nblocks_clinic
         this_block = get_block(blocks_clinic(iblock),iblock)
         do j=this_block%jb,this_block%je
         do i=this_block%ib,this_block%ie
            n = n + 1
            SHR_STREAM_WORK(i,j,iblock) = &
               ndep_sdat%avs(1)%rAttr(ndep_shr_stream_no_ind,n)
         enddo
         enddo
      enddo

      call POP_HaloUpdate(SHR_STREAM_WORK,POP_haloClinic, &
                          POP_gridHorzLocCenter,          &
                          POP_fieldKindScalar, errorCode, &
                          fillValue = 0.0_POP_r8)
      if (errorCode /= POP_Success) then
         call exit_POP(sigAbort, subname /&
            &/ ': error updating halo for Ndep fields')
      endif

      !$OMP PARALLEL DO PRIVATE(iblock)
      do iblock = 1, nblocks_clinic
         where (LAND_MASK(:,:,iblock))
            STF_MODULE(:,:,no3_ind,iblock) = &
               ndep_shr_stream_scale_factor * SHR_STREAM_WORK(:,:,iblock)
         elsewhere
            STF_MODULE(:,:,no3_ind,iblock) = c0
         endwhere
      enddo
      !$OMP END PARALLEL DO

      !
      ! process NH4 flux, store results in SHR_STREAM_WORK array
      ! instead of directly into STF_MODULE
      ! to avoid argument copies in HaloUpdate calls
      !
      n = 0
      do iblock = 1, nblocks_clinic
         this_block = get_block(blocks_clinic(iblock),iblock)
         do j=this_block%jb,this_block%je
         do i=this_block%ib,this_block%ie
            n = n + 1
            SHR_STREAM_WORK(i,j,iblock) = &
               ndep_sdat%avs(1)%rAttr(ndep_shr_stream_nh_ind,n)
         enddo
         enddo
      enddo

      call POP_HaloUpdate(SHR_STREAM_WORK,POP_haloClinic, &
                          POP_gridHorzLocCenter,          &
                          POP_fieldKindScalar, errorCode, &
                          fillValue = 0.0_POP_r8)
      if (errorCode /= POP_Success) then
         call exit_POP(sigAbort, subname /&
            &/ ': error updating halo for Ndep fields')
      endif

      !$OMP PARALLEL DO PRIVATE(iblock)
      do iblock = 1, nblocks_clinic
         where (LAND_MASK(:,:,iblock))
            STF_MODULE(:,:,nh4_ind,iblock) = &
               ndep_shr_stream_scale_factor * SHR_STREAM_WORK(:,:,iblock)
         elsewhere
            STF_MODULE(:,:,nh4_ind,iblock) = c0
         endwhere
      enddo
      !$OMP END PARALLEL DO

   endif
#endif

   call timer_stop(ecosys_sflux_timer)

!-----------------------------------------------------------------------
!EOC

 end subroutine ecosys_set_sflux

!*****************************************************************************
!BOP
! !IROUTINE: SCHMIDT_O2
! !INTERFACE:

 function SCHMIDT_O2(SST, LAND_MASK)

! !DESCRIPTION:
!  Compute Schmidt number of O2 in seawater as function of SST
!  where LAND_MASK is true. Give zero where LAND_MASK is false.
!
!  ref : Keeling et al, Global Biogeochem. Cycles, Vol. 12,
!        No. 1, pp. 141-163, March 1998
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block), intent(in) :: SST

   logical (log_kind), dimension(nx_block,ny_block), intent(in) :: &
      LAND_MASK

! !OUTPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block) :: SCHMIDT_O2

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   real (r8), parameter :: &
      a = 1638.0_r8, &
      b = 81.83_r8, &
      c = 1.483_r8, &
      d = 0.008004_r8

!-----------------------------------------------------------------------

   where (LAND_MASK)
      SCHMIDT_O2 = a + SST * (-b + SST * (c + SST * (-d)))
   elsewhere
      SCHMIDT_O2 = c0
   end where

!-----------------------------------------------------------------------
!EOC

 end function SCHMIDT_O2

!*****************************************************************************
!BOP
! !IROUTINE: O2SAT
! !INTERFACE:

 function O2SAT(SST, SSS, LAND_MASK)

! !DESCRIPTION:
!
!  Computes oxygen saturation concentration at 1 atm total pressure
!  in mmol/m^3 given the temperature (t, in deg C) and the salinity (s,
!  in permil) where LAND_MASK is true. Give zero where LAND_MASK is false.
!
!  FROM GARCIA AND GORDON (1992), LIMNOLOGY and OCEANOGRAPHY.
!  THE FORMULA USED IS FROM PAGE 1310, EQUATION (8).
!
!  *** NOTE: THE "A_3*TS^2" TERM (IN THE PAPER) IS INCORRECT. ***
!  *** IT SHOULD NOT BE THERE.                                ***
!
!  O2SAT IS DEFINED BETWEEN T(freezing) <= T <= 40(deg C) AND
!  0 permil <= S <= 42 permil
!  CHECK VALUE:  T = 10.0 deg C, S = 35.0 permil,
!  O2SAT = 282.015 mmol/m^3
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block), intent(in) :: &
      SST, & ! sea surface temperature (C)
      SSS    ! sea surface salinity (psu)

   logical (log_kind), dimension(nx_block,ny_block), intent(in) :: &
      LAND_MASK

! !OUTPUT PARAMETERS:

    real (r8), dimension(nx_block,ny_block) :: O2SAT

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   real (r8), dimension(nx_block,ny_block) :: TS

!-----------------------------------------------------------------------
!  coefficients in expansion
!-----------------------------------------------------------------------

   real (r8), parameter :: &
      a_0 = 2.00907_r8, &
      a_1 = 3.22014_r8, &
      a_2 = 4.05010_r8, &
      a_3 = 4.94457_r8, &
      a_4 = -2.56847E-1_r8, &
      a_5 = 3.88767_r8, &
      b_0 = -6.24523E-3_r8, &
      b_1 = -7.37614E-3_r8, &
      b_2 = -1.03410E-2_r8, &
      b_3 = -8.17083E-3_r8, &
      c_0 = -4.88682E-7_r8

   where (LAND_MASK)
      TS = log( ((T0_Kelvin+25.0_r8) - SST) / (T0_Kelvin + SST) )

      O2SAT = exp(a_0+TS*(a_1+TS*(a_2+TS*(a_3+TS*(a_4+TS*a_5)))) + &
         SSS*( (b_0+TS*(b_1+TS*(b_2+TS*b_3))) + SSS*c_0 ))
   elsewhere
      O2SAT = c0
   end where

!-----------------------------------------------------------------------
!  Convert from ml/l to mmol/m^3
!-----------------------------------------------------------------------

   O2SAT = O2SAT / 0.0223916_r8

!-----------------------------------------------------------------------
!EOC

 end function O2SAT

!*****************************************************************************
!BOP
! !IROUTINE: SCHMIDT_CO2
! !INTERFACE:

 function SCHMIDT_CO2(SST, LAND_MASK)

! !DESCRIPTION:
!  Compute Schmidt number of CO2 in seawater as function of SST
!  where LAND_MASK is true. Give zero where LAND_MASK is false.
!
!  ref : Wanninkhof, J. Geophys. Res, Vol. 97, No. C5,
!  pp. 7373-7382, May 15, 1992
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block), intent(in) :: SST

   logical (log_kind), dimension(nx_block,ny_block), intent(in) :: &
      LAND_MASK

! !OUTPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block) :: SCHMIDT_CO2

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   real (r8), parameter :: &
      a = 2073.1_r8, &
      b = 125.62_r8, &
      c = 3.6276_r8, &
      d = 0.043219_r8

   where (LAND_MASK)
      SCHMIDT_CO2 = a + SST * (-b + SST * (c + SST * (-d)))
   elsewhere
      SCHMIDT_CO2 = c0
   end where

!-----------------------------------------------------------------------
!EOC

 end function SCHMIDT_CO2

!*****************************************************************************
!BOP
! !IROUTINE: ecosys_tavg_forcing
! !INTERFACE:

 subroutine ecosys_tavg_forcing(STF_MODULE)

! !DESCRIPTION:
!  Accumulate non-standard forcing related tavg variables.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

  real (r8), dimension(nx_block,ny_block,ecosys_tracer_cnt,max_blocks_clinic), &
     intent(in) :: STF_MODULE

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      iblock              ! block loop index

!-----------------------------------------------------------------------
!  multiply IRON, DUST fluxes by mpercm (.01) to convert from model
!    units (cm/s)(mmol/m^3) to mmol/s/m^2
!-----------------------------------------------------------------------

   !$OMP PARALLEL DO PRIVATE(iblock)

   do iblock = 1,nblocks_clinic

      call accumulate_tavg_field(STF_MODULE(:,:,fe_ind,iblock)*mpercm, &
                                 tavg_IRON_FLUX,iblock,1)
      call accumulate_tavg_field(dust_FLUX_IN(:,:,iblock)*mpercm,  &
                                    tavg_DUST_FLUX,iblock,1)
      call accumulate_tavg_field(STF_MODULE(:,:,no3_ind,iblock),  &
                                 tavg_NOx_FLUX,iblock,1)
      call accumulate_tavg_field(STF_MODULE(:,:,nh4_ind,iblock),  &
                                 tavg_NHy_FLUX,iblock,1)

!-----------------------------------------------------------------------
!  now do components of flux calculations saved in hack array
!-----------------------------------------------------------------------

      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,1,iblock),  &
                                 tavg_ECOSYS_IFRAC,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,1,iblock),  &
                                 tavg_ECOSYS_IFRAC_2,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,2,iblock),  &
                                 tavg_ECOSYS_XKW,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,2,iblock),  &
                                 tavg_ECOSYS_XKW_2,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,3,iblock),  &
                                 tavg_ECOSYS_ATM_PRESS,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,4,iblock),  &
                                 tavg_PV_O2,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,5,iblock),  &
                                 tavg_SCHMIDT_O2,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,6,iblock),  &
                                 tavg_O2SAT,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,7,iblock),  &
                                 tavg_CO2STAR,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,8,iblock),  &
                                 tavg_DCO2STAR,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,9,iblock),  &
                                 tavg_pCO2SURF,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,10,iblock),  &
                                 tavg_DpCO2,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,10,iblock),  &
                                 tavg_DpCO2_2,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,11,iblock),  &
                                 tavg_PV_CO2,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,12,iblock),  &
                                 tavg_SCHMIDT_CO2,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,13,iblock),  &
                                 tavg_DIC_GAS_FLUX,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,13,iblock),  &
                                 tavg_DIC_GAS_FLUX_2,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,14,iblock),  &
                                 tavg_PH,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,15,iblock),  &
                                 tavg_ATM_CO2,iblock,1)
      call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,16,iblock),  &
                                 tavg_O2_GAS_FLUX_2,iblock,1)

   end do

   !$OMP END PARALLEL DO

!-----------------------------------------------------------------------
!EOC

 end subroutine ecosys_tavg_forcing

!*****************************************************************************
!BOP
! !IROUTINE: comp_surf_avg
! !INTERFACE:

 subroutine comp_surf_avg(SURF_VALS_OLD,SURF_VALS_CUR)

! !DESCRIPTION:
!  compute average surface tracer values
!
!  avg = sum(SURF_VAL*TAREA) / sum(TAREA)
!  with the sum taken over ocean points only
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

  real (r8), dimension(nx_block,ny_block,ecosys_tracer_cnt,max_blocks_clinic), &
      intent(in) :: SURF_VALS_OLD, SURF_VALS_CUR

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      n,       & ! tracer index
      iblock,  & ! block index
      dcount,  & ! diag counter
      ib,ie,jb,je

   real (r8), dimension(max_blocks_clinic,ecosys_tracer_cnt) :: &
      local_sums ! array for holding block sums of each diagnostic

   real (r8) :: &
      sum_tmp ! temp for local sum

   real (r8), dimension(nx_block,ny_block) :: &
      WORK1, &! local work space
      TFACT   ! factor for normalizing sums

   type (block) :: &
      this_block ! block information for current block

!-----------------------------------------------------------------------

   local_sums = c0

!jw   !$OMP PARALLEL DO PRIVATE(iblock,this_block,ib,ie,jb,je,TFACT,n,WORK1)
   do iblock = 1,nblocks_clinic
      this_block = get_block(blocks_clinic(iblock),iblock)
      ib = this_block%ib
      ie = this_block%ie
      jb = this_block%jb
      je = this_block%je
      TFACT = TAREA(:,:,iblock)*RCALCT(:,:,iblock)

      do n = 1, ecosys_tracer_cnt
         if (vflux_flag(n)) then
            WORK1 = p5*(SURF_VALS_OLD(:,:,n,iblock) + &
                        SURF_VALS_CUR(:,:,n,iblock))*TFACT
            local_sums(iblock,n) = sum(WORK1(ib:ie,jb:je))
         endif
      end do
   end do
!jw   !$OMP END PARALLEL DO

   do n = 1, ecosys_tracer_cnt
      if (vflux_flag(n)) then
         sum_tmp = sum(local_sums(:,n))
         surf_avg(n) = global_sum(sum_tmp,distrb_clinic)/area_t
      endif
   end do

   if(my_task == master_task) then
      write(stdout,*)' Calculating surface tracer averages'
      do n = 1, ecosys_tracer_cnt
         if (vflux_flag(n)) then
            write(stdout,*) n, surf_avg(n)
         endif
      end do
   endif

!-----------------------------------------------------------------------
!EOC

 end subroutine comp_surf_avg

!*****************************************************************************
!BOP
! !IROUTINE: ecosys_write_restart
! !INTERFACE:

 subroutine ecosys_write_restart(restart_file, action)

! !DESCRIPTION:
!  write auxiliary fields & scalars to restart files
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   character(*), intent(in) :: action

! !INPUT/OUTPUT PARAMETERS:

   type (datafile), intent (inout)  :: restart_file

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   character (char_len) :: &
      short_name  ! tracer name temporaries

   type (io_dim) :: &
      i_dim, j_dim ! dimension descriptors

   integer (int_kind) :: n

   type (io_field_desc), save :: PH_SURF

!-----------------------------------------------------------------------

   if (trim(action) == 'add_attrib_file') then
      short_name = char_blank
      do n=1,ecosys_tracer_cnt
         if (vflux_flag(n)) then
            short_name = 'surf_avg_' /&
                      &/ ind_name_table(n)%name
            call add_attrib_file(restart_file,trim(short_name),surf_avg(n))
         endif
      end do
   endif

   if (trim(action) == 'define') then
      i_dim = construct_io_dim('i', nx_global)
      j_dim = construct_io_dim('j', ny_global)

      PH_SURF = construct_io_field('PH_SURF', i_dim, j_dim,     &
                   long_name='surface pH at current time',      &
                   units    ='pH', grid_loc ='2110',            &
                   field_loc = field_loc_center,                &
                   field_type = field_type_scalar,              &
                   d2d_array = PH_PREV)
      call data_set (restart_file, 'define', PH_SURF)
   endif

   if (trim(action) == 'write') then
      call data_set (restart_file, 'write', PH_SURF)
   endif

!-----------------------------------------------------------------------
!EOC

 end subroutine ecosys_write_restart

!*****************************************************************************
!BOP
! !IROUTINE: ecosys_tracer_ref_val
! !INTERFACE:

 function ecosys_tracer_ref_val(ind)

! !DESCRIPTION:
!  return reference value for tracers using virtual fluxes
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   integer (int_kind), intent(in) :: ind

! !OUTPUT PARAMETERS:

   real (r8) :: ecosys_tracer_ref_val

!EOP
!BOC
!-----------------------------------------------------------------------

   if (vflux_flag(ind)) then
      ecosys_tracer_ref_val = surf_avg(ind)
   else
      ecosys_tracer_ref_val = c0
   endif

!-----------------------------------------------------------------------
!EOC

 end function ecosys_tracer_ref_val

!***********************************************************************

 end module ecosys_mod

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
